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A Fokker-Planck Study of Dense Rotating Stelleir Clusters 

Abstract 

The dynamical evolution of dense stellar systems is simulated using a two-dimensional 
Fokker-Planck method, with the goal of providing a model for the formation of supermassive 
stars which could serve as seed objects for the supermassive black holes of quasars. This 
work follows and expands on earlier one-dimensional studies of spherical clusters of main- 
sequence stars. The two-dimensional approach allows for the study of rotating systems, 
as would be expected due to cosmological tidal torquing; other physical effects included 
are collisional mergers of individual stars and a bulk stellar bar perturbation in the sys- 
tem's gravitational potential. The 3 Myr main-sequence lifetime for large stars provides 
an upper limit on the allowed simulation times. Two general classes of initial systems are 
studied: Plummer spheres, which represent stellar clusters, and "7 = 0" spheres, which 
model galactic spheroids. 

At the initial densities of the modeled systems, mass segregation and runaway stellar 
collisions alone are insufficient to induce core collapse within the main-sequence lifetime 
limit, if no bar perturbation is included. However, core collapse is not a requirement for 
the formation of a massive object: the choice of stellar initial mass function (IMF) is 
found to play a crucial role. When using an IMF similar to that observed for dense stellar 
clusters (weighted towards high masses but with a high-mass cutoff of M^^x ^ ISOMq) the 
simulations presented here show, in all cases, that the stellar system forms massive (25OM0) 
objects by collisional mergers of lower-mass stars; in almost all such cases the presence of a 
stellar bar allows for sufficient additional outward transport of angular momentum that a 
core-collapse state is reached with corresponding further increase in the rate of formation 
of massive objects. In contrast, simulations using an IMF similar to that observed for field 
stars in general (which is weighted more towards lower masses) produce no massive objects, 
and reach core collapse only for initial models which represent the highest-density galactic 
spheriods. 

Possible extensions of the work presented here include continuing to track stellar 
populations after they evolve off the main sequence, and allowing for a (possibly changing) 
nonspherical component to the overall system potential. 
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values into the calculation produced identical curves. Similarly, whether the 
VLj were found by integrating over the orbit (using (j2.2p or (j2.3p ) or by taking 
Q.j = ^ (|2.5p made no difference 

4.4 Qi (upper plot) and Q2 vs. / for an isochrone potential, at a constant energy. 
Similar to Fig. 14.31 but here the numerically-found curves show both the 
semi-analytic- and fully-numeric-potential cases 

4.5 Orbital angle wi{r) interpolated from values calculated using (12.61) v. the 
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6 = 71/2 for the most-curved 
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model. Full results for these simulations are given in Chapter [5l 

4.12 Comparison of different flnite-differencing and boundary-condition schemes. 
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analytic solution is identical to the t = Q curve at all times). Here f{x,y) is 
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4.14 Comparison of results of using different-size (/, J) grids. The break in the 
otherwise smooth 34x34 curve is a spurious artifact of the output-calculating 
procedure and does not affect the model's further evolution. All other pa- 
rameters are the same for each run, and each is constrained to having the 
same grid as the others at each timestep. (These runs used grids with 
50 points, a maximum I value of 9 and a 3-component Miller-Scalo-style IMF.) 

4.15 Comparison of runs with grids containing 44, 50 and 60 points; other 
parameters (timestep, action-space gridsize etc.) were identical for each of 
the 3 cases. In order to do a direct comparsion, each run was prevented 
from adjusting its grid dynamically. These runs were for model E4B with 
^max = 4 and an initial rotation parameter A = 0.05. The mass spectrum was 
90% IMq and 10% 2Mq 

4.16 Comparison of runs with grids containing 40, 50, 60 and 74 points using 
a 40x40 (I, J) grid, /max = 3 and identical timesteps in each case. Dynamic 
7-2 grids were enabled, but were similar in all 3 cases. These runs were for a 
7-mass component cluster (mass range 1 — 32Mq in 7 bins) with a moderate 
amount of rotation (Aq = 0.015), Kroupa-style IMF and with 1% of the 
cluster mass in a stellar bar 

4.17 Comparison of the use of different-size timesteps. Runs are with r^-grids of 
50 points; the At = 2.4 case continued until t ~ 73, as can be seen in Fig. 
14.151 Initial central relaxation time was trq{0) = 25.6 Myr 

4.18 Comparison of the use of different-size timesteps. Runs are of model E4B 
with an r^-grid of 74 points, a 40x40 (I, J) grid, a maximum I of 2 and an 
initial rotation given by Aq = 0.05 and a Kroupa IMF with 10 stellar mass 
bins in the range 1 — 125M0. Initial central relaxation time for the IMq stars 
that comprised 51% of the cluster was iri(O) = 158 Myr; for comparison the 
25M0 stars' mass fraction was 2.4% with tr7(0) = 0.23 Myr 

4.19 Similar to Fig. 14.181 but showing the full range of the At = 1.52 run. . . . 

4.20 Demonstration of the use of an /max value of 3, 4, 5 and 6 in the diffusion 
coefficient expansion ()3.23p (and by extension in the expansion of the poten- 
tial (j3.14p ): this corresponds to a total number of expansion terms of 4, 6, 
9 and 12 respectively. This plot is of a rotating 2-component cluster similar 
to those of Figures 14.151 and 14.171 Each run was constrained to use the same 
timestep and the same r^-grid as the others, to allow for direct comparison. 

4.21 Comparison of different numbers of expansion terms and gridpoints for 
the case including a bar perturbation. The 50-point, quadrupole-only run 
continued on a fairly linear path (not shown) until it reached p ~ 2 x 10^ 
at t = 0.2. This run is for a 2-component cluster with 10% 2Mq and 90% 
IMq stars, 4% of which comprised the bar. The rotational parameter was 
A = 0.075. A run with /max = 4 {i.e., 6 expansion terms) was nearly identical 
to the one shown using /max = 3, although it ended earlier 

4.22 Comparison of different numbers of expansion terms for a "production" run: 
model E4B with a Kroupa IMF, initial rotation parameter Aq = 0.05 and a 
stellar bar mass fraction of 1%. The stellar mass range is 1 — 125 Mq. . . 
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4.23 Comparison of different numbers of expansion terms for a "production" run: 
model E2A with a Kroupa IMF, initial rotation parameter Aq = 0.05 and a 
stellar bar mass fraction of 1%. The stellar mass range is 1 — 125 Mq. . . 

4.24 Comparison of central density for runs in which the stellar mass spectrum 
is split into the number of geometrically-spaced mass bins shown. All are 
model E4B with 50 points in the grid, a 34x34 (/, J) grid, and an overall 
mass range of 1 — 32M0. The IMF is Salpeter-like, but with the uppermost 
mass bin artifically overpopulated in order to accentuate any differences in 
the cases 

4.25 Similar to Fig. I4.24l but showing the central density of the uppermost {mq = 
32Mq) stellar mass bin only 

4.26 Comparison of central density for runs in which the stellar mass spectrum 
is split into the number of geometrically-spaced mass bins shown. All are 
model E4B with 50 points in the grid, a 34x34 (/, J) grid, and an overall 
mass range of 1 — 32M0. The IMF is Kroupa-like, similar to that used in 
the later full simulations. Rotation and a stellar bar are included, with a bar 
mass fraction of 1% and an initial rotation parameter of Aq = 0.016 

4.27 Similar to Fig. 14.261 but with stellar mergers included in the simulation. 
(The 9-bin case shows slightly lower central density than the others due to 
it following a somewhat different series of values for the innermost point 
of its dynamic grid; this was unusual and did not occur during the full 
simulations presented in Chapter [5] except where noted.) 

4.28 Ratio of maximum possible amount of binary heating Kb to overall (grav- 
itational + kinetic) system energy \Etot\ for a variety of models. All used 
the Arches-style IMF except for model "E2A/Kroupa" , and all had an ini- 
tial Plummer-sphere distribution except for "G2A" , which started as a 7 = 
sphere. A logarithmic scale is used so that different models' timescales can be 
plotted together; exactly flat lines at the start of each model's plot indicate 
the first timestep in the calculation, before which binary heating is zero. 

5.1 Central density v. time for model E4B using the Kroupa IMF with a stellar 
mass range of 1 — 125M0 and /max = 2. Note that the nonrotating (Aq = 0) 
case has an initial p{0) slightly higher than the cases that include rotation; 
this is a product of the method for introducing rotation into the initial dis- 
tribution function of the system 

5.2 Similar to Fig. I5.1l but showing only the highest-mass stars (mg = I25M0). 
The upturn at the end of the "with mergers" curve is for the final timestep 
only and is a result of a numerical instability which casued the simulation to 
end early 
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5.3 Density p{r) v. radial distance r from the cluster center for the start, mid- 
point and end of the simulation for model E4B with a Kroupa IMF. Stellar 
mergers were enabled, as was a 10% stellar bar. Rotation parameter A = 0.05 
and /max = 2. (This corresponds to the uppermost plot of Fig. 15.11 ) Note 
that the decrease in p{r) with time for large r is only exhibited when a stellar 

bar is present [t^ 

5.4 Central density v. time for model E4B using the Arches IMF with a stellar 
mass range of 2 — 125Mq and /max = 2 [s^ 

5.5 Similar to Fig. 15.41 but showing only the highest-mass stars (mg = 125Mq). [s^ 

5.6 Density p{r) v. radial distance r from the cluster center for the start, mid- 
point and end of the simulation for model E4B with an Arches-style IMF. 
Stellar mergers were enabled, as was a 10% stellar bar; rotation parameter 

A = 0.049 and /max = 2. (This corresponds to the leftmost plot of Fig. 15.41 ) [81 

5.7 Density p{r) v. radial distance r from the cluster center for the start and end 
of the simulation for model E4B/ Arches, as well as at two similar times as 
those shown in Fig. 15.61 for comparison. Neither stellar mergers nor a stellar 
bar were enabled (in contrast with Fig. 15. 6p . Rotation parameter A = 0.049 

and /max = 2 [Sl 

5.8 Central velocity dispersion (Ti(0) v. time of the lowest-mass {2Mq) stars 
in model E4B with an Arches IMF, for nonrotating, rotating and with-bar 
cases. (The case shown with a stellar bar also had stellar mergers enabled, 

but a run without mergers was similar.) [s^ 

5.9 Median central relaxation time tr{0) v. time of the lowest-mass (IMq) stars 
in model E4B with an Arches IMF, for nonrotating, rotating and with-bar 
cases. (The case shown with a stellar bar also had stellar mergers enabled, 

but a run without mergers was similar.) [s^ 

5.10 Central density v. time for model E2A with Arches-style IMF, /max = 3 and 

a stellar mass range of 2 — 125Mq [84 

5.11 Central density v. time for model E2B with /max = 2 and using the Arches 
IMF with a stellar mass range of 2-125M0. Note that the main-sequence 
lifetime of the most massive stars is 3 Myr and so the region of the graph 
beyond t = 3 is nonphysical [84 

5.12 Detail of the first 3 Myr of Fig. ISTTl [85 

5.13 Similar to Fig. 15.111 but including /max = 3 and an 8% bar case. All have 
A = 0.049 and do not include stellar mergers. (Simulations performed with 
mergers enabled yielded almost-identical runs of p{0) to those shown in this 
figure.) Note that the "5% bar" plot here is the same as that shown in Fig. 
15.111 The gap in one plot is due to a glitch in the output routine, as described 

in Chapter H [sE 
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5.14 Density p(r) v. radial distance r from the cluster center for the start, mid- 
point and end of the simulation for model E2B with an Arches-style IMF. 
Stellar mergers were enabled, as was an 8% stellar bar; rotation parameter 
A = 0.049, and /max = 3. This corresponds to the leftmost plot of Fig. 15.131 
The decrease in p{r) with time for large r is only seen when a stellar bar is 
present 

5.15 Effect of artificially increasing the collisional merger rate. Plotted is model 
E4B using a Kroupa IMF with /max = 2 and a 10% stellar bar. Also see 
Table EH 

5.16 Central velocity dispersion v. time for the lowest-mass (IMq) stars in the 
7 = sphere models with a Kroupa IMF. All have A ~ 0.05. Plots of runs 
with and without stellar mergers were identical to each other; differences 
within each model are due to the presence or absence of a stellar bar. . . . 

5.17 Central relaxation time for the lowest-mass {IMq) stars in the 7 = sphere 
models with a Kroupa IMF. All have A ~ 0.05 

5.18 Central density v. time for 7 = sphere model G2A with a Kroupa IMF, 
/max = 4 and A = 0.051. The nonrotating (A = 0) case gave similar results 
but was less stable and not plotted 

5.19 Density p(r) v. radial distance r from the cluster center for the start, mid- 
point and end of the simulation for model G2A with a Kroupa IMF. Neither 
rotation nor stellar mergers were enabled, and /max = 4; this corresponds to 
the solid-line plot of Fig. 15.181 The behavior of p{0) v. r for the other two 
cases from Fig. 15.181 is almost identical to that plotted here and are so not 
shown 

5.20 Central density v. time for 7 = sphere model G3A with a Kroupa IMF, 
/max = 3, A = 0.046 and a 7.7% stellar bar. The gap in the no-merger plot 
is due to a glitch in the output routine. Runs without a stellar bar did not 
yield stable simulations and are not plotted 

5.21 Density p{r) v. radial distance r from the cluster center for the start, 
midpoint and end of the simulation for model G3A with a Kroupa IMF, 
A = 0.046, a 7.7% stellar bar, and stellar mergers enabled. This corresponds 
to the dashed plot of Fig. KM 

5.22 Central density v. time for 7 = sphere model G3C with a Kroupa IMF, 
/max = 2, and A = 0.046. The nonrotating (A = 0) case did not yield a stable 
simulation and is not plotted 

5.23 Density p(r) v. radial distance r from the cluster center for the start, 
midpoint and end of the simulation for model G3C with a Kroupa IMF, 
A = 0.046, a 7.7% stellar bar, and stellar mergers enabled. This corresponds 
to the short-dashed plot of Fig. 15.221 

5.24 Central density v. time for 7 = sphere model G2A with an Arches-style 
IMF, /max = 3 and A = 0.051. Two of the runs continued beyond t = 0.07 
Myr but became unstable; those portions are not shown 
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2.1 The set of initial condition models. For purposes of the numerical simulation, 
the quantities taken as fundamental parameters are total mass M and core 
radius rcore- The central density /5(0), velocity dispersion o"o(0) and relaxation 
time tr'(O) listed here are the resulting analytic values for a Plummer sphere 
distribution. Model names are taken from Quinlan and Shapiro [l] 

5.1 Models which yielded stable results, along with possible analogous astronom- 
ical systems 

5.2 List of models which yielded stable simulations using the Kroupa IMF, with 
central-density and merger-rate results shown at various times ti in each 
simulation (usually the endpoint, or chosen for comparison with another). 
Times are in Myr, central densities in Mq/pc'^, and each simulation had a 
rotation parameter A ~ 0.05. 6*250 is the calculated rate at which 250Mq 
stars would be expected to be produced via collisional mergers of 125M0 
stars, in Myr^-*^. S is the factor by which the central density of the most- 
massive stars increased as a result of mass segregation, relative to the overall 
increase in central density; S' is the same factor but also including the effects 
of stellar mergers. For the G2A and G3C models, the presence or absence of 
a stellar bar had little effect on the system's evolution. (Model G3A was not 
run without a stellar bar but it was also expected to have little effect.) Note 
that the central densities listed here are for the initial models with rotation, 
and so are lower than those given in Table [2Tl Starting time to = 0. ... 

5.3 Similar to Table 15.21 but for simulations using an Arches-style IMF. [*] in- 
dicates that the value of ^250(^1) calculated for model G2A without a bar is 
considered to be numerically suspect 

5.4 Addendum to Table 15.31 showing results of simulations allowed to proceed 
beyond the 3 Myr lifetime of the largest main-sequence stars 

5.5 Effect of artificially increasing the collisional merger rate by a factor of 5 in 
order to account for the effect of tidal-capture binaries, as described in ^5.2.61 
Listed are results for model E4B with Zmax = 2 and a 10% stellar bar. The 
factor = 1 entry has slight differences from the values shown in Table 15.31 
due to different timestep sizes being used for the simulations listed here. . . 

B.l List of frequently-used subscripts 
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B.2 Symbols first introduced in Chapter [2l listed by section in which each was 
first introduced. These denote properties of individual stars in the system, 
not bulk properties. Bulk properties based on summing over the individual 
stars' values are denoted by subscripts, for example total energy of the stellar 
system E'tot or the rotational contribution to the total angular momentum 
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section are included lid 



XVll 



Chapter 1 

Introduction 



"The challenge of simulating a dense star cluster with a million stars is formidable, 
because of the enormous ranges in spatial and temporal scales that have to be 
modeled simultaneously." [2] 

1.1 Astrophysical Motivation: Quasar Massive Black Holes 

The observation of high-redshift quasars, i.e., those with redshift z > 6, implies a problem 
for understanding how they formed. The prevailing model for the underlying astronomical 
object is that of a supermassive black hole accreting material gravitationally [3]. Each 
high-redshift quasar thus must have had a supermassive black hole in place within the first 
10^ years after the big bang, which turns out to be a rather strong physical constraint. 
Observations of high-redshift gamma-ray bursts and massive galaxies place similar limits 
on the formation timescales of the first massive stars [4J and of the first massive galaxies 

Of the two possible basic building blocks for a massive object in the early universe - gas, and 
stars that have already coalesced from the gas - direct formation from the collapse of a region 
of gas may seem to have certain advantages, in that it avoids the intermediate step of forming 
stars. But the required gas cooling times are barely consistent with current structure 
formation models (and even then require reliance on "extremely rare" overdensities) and 
additionally the centrifugal barrier that arises due to tidally-induced rotation is a "key 
obstacle" which must be overcome by any model of massive black hole formation [6j. As 
will be described in more detail below, few previous studies have examined the evolution 
and possible core collapse of rotating dense stellar clusters, and none have modeled whether 
such systems can form massive stars which could then serve as seeds for quasar black holes. 

The task of forming a massive object from stellar- mass objects faces its own hurdles: even 
in dense environments stars are less affected by collisions than is gas (although stellar 
collisions have the benefit of not necessarily causing shock heating) and so energy and 
angular momentum transfer due to collisional processes are not as effective a prioi for stars 
as they are for gas. However, stellar systems can overcome this handicap through the 
formation of a large, coherent dynamical subsystem - the most observationally prevalent of 
which is a "stellar bar" - which can act as a bulk perturbation in the gravitational potential 
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of the system and so can effectively transport angular momentum outwards, thus allowing 
the system's core to contract further. This mechanism is described further in §1.1.31 and 
developed in detail in ^3.51 

The goal of the current project is to generalize the approach of prior researchers, who de- 
veloped one-dimensional models which simulated nonrotating stellar systems {e.g., Quinlan 
and Shapiro [7] and other references given below), to two dimensions in order to address 
the questions: What about rotation, does it inhibit collapse and can it be overcome? The 
anticipated answer is that, despite initial rotational support against collapse, the presence of 
a stellar bar will allow sufficient transport angular momentum outwards so that collisional 
stellar mergers and mass segregation can produce a massive object (perhaps ~ 10^^'^ Mq) 
in the core of a dense stellar cluster such as could be found at the center of a newly-forming 
galaxy. Such an object could then undergo growth via accretion to reach supermassive 
size (~ IO^^^Mq) within a Hubble time [8j and be observed as a quasar or active galactic 
nucleus, thus bringing together the formation of early-universe stars, galaxies, and massive 
collapsed objects. 

The remainder of this chapter expands on each of the above ideas, starting with dense stellar 
systems, and moving on to how massive black holes may form in such stellar clusters. How 
rotation impacts the system is then discussed. Finally an overview of the motivation for 
and development of the numerical simulation technique used is given. 

1.1.1 Stellar Systems 

A primary motivation for this study was the work of Quinlan and Shapiro [l], who developed 
a Fokker-Planck model of an idealized nonrotating, spherically symmetric, dense cluster of 
compact stars, and found "rapid buildup of massive black holes in the cluster core resulting 
from successive binary mergers and mass segregation." Subsequently, Quinlan and Shapiro 
studied clusters of solar-mass main sequence stars and found that it was "remarkably easy 
for massive stars to form through multiple stellar mergers in dense galactic nuclei" [7]. 
After performing a complemetary line of simulations using the specialized GRAPE N-body 
computing cluster, Portegies Zwart and McMillan [9j speculated that a sufficiently-dense 
IO^Mq stellar cluster which initially formed ~ 30 pc from the Galactic center can spiral 
inward due to dynamical friction before being disrupted by the Milky Way's tidal field. 
While a less-dense cluster would be disrupted before reaching the Galactic center and so 
merely contribute its stars to the Galactic bulge, the more-dense cluster dissolves closer to 
the center of the Galaxy and leaves behind its central massive black hole which can then 
grow to supermassive size by accretion. 

In early theoretical work on stellar systems Begelman and Rees [10] found that "runaway 
coalescence" of stars (i.e., stellar mergers combined with mass segregation) could lead to 
the formation of a central massive object. Again, the system's evolution was a race between 
coalescence and disruption, the latter in this case being due to gas released by earlier 
stellar evolution and collisions. Lower-density clusters were susceptible to disruption and 
dissipation, while more-dense systems retained most of the gas which went to form new stars. 
Stellar clusters were thus explicitly linked to active galactic nuclei: "Dense star clusters may 
be responsible for some of the low-level manifestations of activity in galatic nuclei; but they 
are probably merely precursor stages of the more spectacular quasar-type phenomena, which 
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develop after a massive object has formed. The authors note that runaway stellar merging 
combined with mass-segregation provides "one of the quickest routes to the formation of a 
massive object in a dense stellar system". 

Star Cluster Observations 

Dense stellar clusters are fairly common in the nuclei of large galaxies: 75% or more of 
late-type (classification Scd-Sm) spiral galaxies have nuclear clusters, as do at least 50% of 
earlier type (Sa-Sc). These nuclear clusters are as compact as globular clusters found in the 
Milky Way, with a typical (half-light) radius of 2 — 5 pc (Boker [11], and references therein), 
but they are massive. The dynamical mass of a nuclear cluster is in the range IO^-IO^Mq, 
i.e., at the very high end of the globular cluster mass function - and similar to the cluster 
masses used here. Boker concludes that nuclear clusters are "an intriguing environment for 
the formation of massive black holes because of their extreme stellar density" . 

An impediment to building massive black holes in stellar clusters is that the initial mass 
function (IMF) which describes the initial distribution of stars across the possible range of 
masses is observed to have a cutoff at around 15OM0 [12], as described in more detail in 
§2.6.31 Thus there is a gap between the IMF's upper cutoff and the 25OM0 minimum stellar 
mass required to leave behind a remnant massive black hole after going supernova [13]. If 
black hole formation is to be feasible this gap must be overcome by dynamical processes in 
the cluster, some candidates for which are considered in the next section. 

1.1.2 Scenarios for Massive Black Hole Formation in Dense Clusters 

Observationally, Ebisuzaki et al. [H] note the existence of what they refer to as the "missing 
link" between stars and supermassive black holes: an IMBH seen in or near a compact stellar 
cluster in the center of galaxy M82. They take this as evidence for supermassive black holes 
being created through mergers of IMBHs, which themselves form via the merging of massive 
stars in compact clusters. Coleman Miller and Colbert [15] also note that IMBHs exist in the 
current universe, and that inner galactic bars can contribute to the growth of a < lOOOM© 
black hole; they also postulate SMBHs from early IMBH mergers, perhaps as stellar clusters 
percolate to the galaxy center and merge. 

Approaching the problem from the other direction, taking the existence of high-redshift 
{e.g., z = 6.41) quasars as a given and examining timescales by which a supermassive 
black hole (SMBH) may grow from an earlier IMBH, Tyler et al. [16] conclude that a 
supermassive star (whether formed as such or as a stellar-merger product) was a more 
likely SMBH progenitor than any of a primordial IMBH, a large Population III star, or 
the merger product of smaller Population-Ill black hole remants. Each of latter would be 
required to form very early and at the extreme upper end of their physically possible mass 
ranges, followed by continuous Eddington-limit accretion. Again, they did not consider 
rotation. 

The overall difficulty of explaining high-redshift quasar black holes is summarized by Diichting 
|17j : "The existence of sufficient primordial black holes would unaccountably close the uni- 
verse" unless the spectrum of masses is fine-tuned. It is "quite difficult to construct viable 

^This citation also includes Rees' prototypical supermassive black hole formation-scenarios flowchart. 
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physical models." Remnants of Population III stars, or IMBHs formed from the collapse of 
early star clusters, also have "serious problems" accounting for z = 6 quasars. 

In reviewing possible formation processes for massive central black holes in nonrotating 
dense stellar clusters with a few thousand to a few million stars, Rasio et al. [18\ found 
two possible paths. If core collapse of the cluster proceeds quickly enough - occurring 
before the typical lifetime of a large main-sequence star, ~ 3 Myr [E]- stellar collisions 
can produce a single supermassive star which then evolves into a intermediate-mass black 
hole (IMBH) of mass 10^ — IO^Mq. Alternatively, if large individual stars go supernova 
before core collapse and leave stellar-mass black hole remnants, the supernova mass loss may 
reverse core collapse somewhat for smaller clusters, while larger systems may still collapse 
relativistically In either case the initial cluster evolution is dominated by the most 
massive stars, i.e., by the high end of the stellar initial mass fuction (IMF). The former 
case of forming a massive object from collisions of main-sequence stars is the object of this 
dissertation. 

While also plausible, modeling the evolution of a cluster of relativistic stellar remnants 
is a sufficiently different as to require a separate, later study. Since main-sequence stellar 
collisions and cluster dynamics "depend crucially on each other" [20] it is unclear what initial 
conditions to assume for a post-main-sequence cluster. And restricting the simulations to 
a main-sequence stellar population also avoids the complications of how to model stellar 
mass loss due to post-main-sequence stellar evolution in which the liberated gas could do 
any of: remaining as gas in the cluster; evaporate from the system; or coalesce into new 
stars, all while the initial stars are transforming into stellar-mass remnants [7J. While also 
interesting, following this stage of the cluster evolution would greatly increase the complexity 
of an already quite complicated model. 

1.1.3 Rotation due to Tidal Torquing 

Rotation is ubiquitous in the universe. A protosystem's expected rotation from being spun 
up by tidal torques due to cosmological perturbations is "almost independent of the per- 
turbation spectrum" and only weakly correlated with the individual overdensity being con- 
sidered pT]. An initial angular momentum is an "important piece of realism ... while 
the probability that the cloud from which a star cluster originates, has zero total angular 
momentum is very small, practically all models have assumed that." [22j. 

Recently, some work has been done on modeling rotating stellar clusters, and the picture 
that emerges is not always intuitive. Rotating isotropic systems were found to be dynami- 
cally stable in the N-body studies performed by Meza [23], but the systems were found to 
be able to rotate "very rapidly" without becoming oblate. Rotation can accelerate the time 
to core collapse in single-component globular clusters [2^ and was found in 2-component 
globular clusters to increase the speed of dynamical evolution via the enhanced outward 
transfer of angular momentum, but only if dynamical friction does not dominate. (How- 
ever, this also results in rotation being a check on mass segregation: massive stars' angular 
momentum J prevents them from sinking as far in as they would otherwise, due to the 
"gravo-gyro instability" |24j.) 

In general, rotation in dense stellar systems is important for galactic nuclei, although the 
vast bulk of previous studies have focused on globular clusters [25] . Most prior work on 
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larger systems has either been to model a specific system and/or to study systems with 
a pre-existing massive black hole {e.g., the modeling of dwarf elliptical galaxy M32 by 
Arabadjis [26j, and references therein). On galactic scales Barnes and Efstathiou [21j did 
find that rotation can inhibit collapse in N-body simulations, but suggest caution with 
analytic predictions of angular momentum evolution in complex stellar systems. 

As rotation does not necessarily induce flattening of the system, the simulations performed 
here assume sphericity for simplicity's sake, and the potential is constructed so that the 
model is self-consistent with that assumption [27]. 

Seeding Black Holes from Low Angular Momentum Material 

One avenue for overcoming the rotational-support problem is to postulate that only the 
least-rotating material contributes to the formation of a seed black hole. Koushiappas et 
al. [6j modeled SMBH seeds from the gas with the least angular momentum in "rare-peak" 
haloes of the early protogalaxies, under the assumption that the distribution of the gas 
has a significant low-angular-momentum tail, which can be extrapolated to lower masses 
and smaller radii from the distribution found for dark matter in cosmological simulations. 
With this assumption they were able to produce seed IMBHs of 10^ Mq for any halo larger 
than 7 x lO^M© - but no black hole at all for smaller haloes, including the bulges of most 
disk galaxies. While giving promising results for large galaxies, they acknowledge that the 
rotational barrier is a key concept and that their model still requires another undetermined 
mechanism for outward transport of angular momentum in order to achieve the gas's collapse 
to a black hole, one which cannot operate for lower-mass haloes and so does not account for 
any continuation of the galaxy/central-object scaling relation to lower-mass galaxies [llj. 
Thus another solution, which does not depend on the specific requirement of unusually rare 
low-rotation initial conditions, is warranted to account for at least some of the observed 
massive black holes in disk galaxies. 

Stellar Bars as Bulk Transporters of Angular Momentum 

Although some angular momentum J is expected to be transported outwards in a rotat- 
ing stellar system via the effects of shear between the faster-rotating inner regions and the 
slower-rotating outer regions, a bulk perturbation in the potential can be much more effec- 
tive than individual stars at transporting angular momentum. In general, any nonaxially- 
symmetric component of the potential in a rotationally supported gravitational system 
fosters J transport over large distances, due to gravitational torques [28J. The amount of 
transport can be as much as an order of magnitude greater than that in a purely spherical 
stellar system [29] . Specifically, a bar-like perturbation in the potential can provide a much 
more effective angular momentum transfer mechanism compared to what independently- 
orbiting stars can effect [30j, with J being transported from inner to outer regions [31j and 
eventually to the stellar system's halo region - the presence of which does not stabilize the 
system against the formation of a bar perturbation in the first place [32]. The bar excites a 
gravitational wake in the stellar distribution: as the bar passes through a region, other stars 
experience a shear force along its orbital path. This has a net drag effect specifically on 
individual stellar orbits in or near resonance with the bar's rotational frequency [33]. This 
"resonance drag" in finite rotating systems is the analogue of classical dynamical friction in 
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an infinite, homogeneous medium |34j . 

Bar perturbations were first observed within the spiral structure of galaxies and are still 
commonly associated with galactic scales. However, inner bars of scale < 100 pc are common 
and form before larger-scale bars [35] . While galactic bars may be largely gas, inner bars are 
mostly stellar and not good at gas transport ^36j. There are dynamical arguments for why 
having both outer and inner stellar bars should be the standard model, with the inner bar 
decoupling from the outer system [S^- Observationally, for example, the SO galaxy NCG 
4621 has a 60 pc counterrotating {i.e., decoupled) core, the "smallest to date" [38]. Thus, 
subject to tests confirming that a given system satisfies the dynamical instability criterion 
for a bar to form (as given in ^3.6p . a stellar bar is invoked here on the scale of the stellar 
cluster being modeled in order to foster angular momentum transport outward from the 
cluster center. 

1.2 Motivation for Technique 

1.2.1 Comparison of Fokker-Planck and N-body Methods 

The two main techniques for modeling collisional stellar systems such as globular clusters, 
open clusters, or galactic nuclei, are through N-body simulations in which particles rep- 
resenting a number of stars are allowed to interact via relevant gravititational and other 
physical processes, and through integration of the Fokker-Planck equation which represents 
the stellar population statistically as a distribution function /. When compared against 
each other on simple (spherical, nonrotating) systems the two techniques agree to a large 
extent ([22]. |39j ) and so the choice of which method to use depends on the requirements of 
a given study. 

N-body simulations have the advantage of being "direct" in that individual particles are 
tracked and can be followed through the entire simulated evolution of the system over time, 
and the physical effects that dominate their behavior can be analyzed. A disadvantage is 
that it is not always apparent how to incorporate a given physical effect into the N-body 
construct in the first place. Nor do individual particles in the simulation correspond to 
individual objects in astronomy {i.e., stars). The typical in current N-body studies is still 
low compared to the number of stars in real astronomical systems, and how to extrapolate 
between "particles" and "stars" is not apparent [22]. Also N-body methods are noisy in the 
sense that a given simulation cannot be relied on to represent the evolution of a system; 
instead one must run a suite of simulations in order to get a statistically robust picture, 
which will be the case as long as N is not up to realistic particle numbers, preventing truly 
star- by-star modeling |25j . 

On the other hand, the statistical nature of the Fokker-Planck equation "captures the 
essence" of continuum mechanics |40j but does not afford direct analysis of how different 
physical processes contribute to the observed behavior of the system. All one knows is how 
the stellar distribution function / evolves over time. However, Fokker-Planck does allow for 
different aspects of the physics to be added or removed at will; thus comparing integrations 
with and without a given effect included can accomplish the same result. This ability to 
include new dynamical aspects in an ad hoc manner permits us to study a discrete spectrum 
of stellar masses {e.g., as previously done by Amaro-Seoane [41]), collisional mergers of stars 
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{e.g., Quinlan and Shapiro [7]) and a bulk perturbation in the gravitational potential such 
as the stellar bar. 

1.2.2 Orbit Averaging and the Third Integral 

As developed in Chapter [3l the Fokker-Planck approximatiorH is appropriate when small 
angle, 2-body scattering dominates the global evolution of the system. When that is the 
case, energy transport can be treated as analogous to heat conduction in a collisional gas 
- producing a "fairly good description of what happens in N-body simulations" ^22j. The 
full Fokker-Planck equation tracks the time-evolution of the system in all 6 dimensions 
of position and velocity (or more generally, momentum) which is both computationally 
prohibitive and usually not necessary for describing the system's evolved properties. 

In order to simplify the problem one can attempt to eliminate the positional dependence 
of the Fokker-Planck equation, leaving only the momentum dependence to describe the 
dynamics. A typical method used to achieve this is oihit averaging in which the Fokker- 
Planck equation and all dynamical quantities on which it depends are, as the name implies, 
averaged over a full stellar orbit before being applied. These orbit-averaged Fokker-Planck 
models are found to "treat very well the diffusion of orbits according to the changes of 
their constants of motion, taking into account the potential and the orbital structure of the 
system in a self-consistent way" |41j . 

One is thus left with a description of a stellar system based on three constants of motion, 
corresponding to the three dimensions of momentum. This is consistent with the strong 
Jeans theorem, which requires in the general case that there be three integrals of the motion, 
even if only two are well determined [39j- Starting with Goodman's 1993 dissertation 
|43j standard practice has been to use energy E and angular momentum component 
as the constants of motion, and to simply ignore any third integral. Comparison to N- 
body results has shown this to be reasonable if flattening of the system is not extreme 
|22) . For such systems, however, it is in principle possible to use as a proxy for the 
third integral, although this would be "extremely difficult numerically and physically" |39j . 
Strictly speaking, ignoring the third integral is theoretically valid if and only if the system's 
gravational potential is spherically symmetric, even if its velocity distribution is not 02] - 
hence the restriction to studying rotating, but gravitationally-spherical systems. 

1.2.3 Choice of Canonical Variables 

Depending on the choice for constants of the motion, the two-dimensional orbit-averaged 
Fokker-Planck equation can be used to study different dynamical apsects of a stellar cluster. 
With the assumption of isotropy of the velocity dispersion, the traditional choice of {E^ J^) 
can represent an axisymmetric rotating system, e.g., Kim et ai.'s single-mass [39] and two- 
mass clusters [23]. Alternatively {E, J^) can be used to study the somewhat simpler problem 
of anisotropy in spherically symmetric systems, as first done by Cohn [13] and more recently 
by Takahashi (L45j, [46j) and Fiestas et al. [25]. The latter problem is simpler in that the 

^Approximation here refers to the idea that the Fokker-Planck equation can be derived by expanding the 
coUisional Boltzmann equation under the assumption of weak encounters, i.e., Sv <^ v, and truncating the 
resulting series after the second-order terms [42) . 
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distribution can be averaged over angular momentum J before coefficients of the Fokker- 
Planck equation are calculated, providing a great savings in complexity and calculation 
time. 

However, once incorporation of ad hoc terms into the Fokker-Planck equation is desired, e.g., 
to represent collisional stellar mergers (as originated by Quinlan and Shapiro [7] in a one- 
dimensional study) energy E is no longer necessarily an integral of the motion. Any change 
in the potential will cause changes in the stellar distribution function when it is expressed 
as f{E,j'^) or as f{E,Jz). What is required are the radial action / and the tangential 
action {i.e., the overall angular momentum) J, which are still adiabatic invariants as the 
potential evolves, and so also is /(/, J). Binney and Tremaine [32] pointed out this exact 
advantage of using the canonical actions as coordinates for direct numerical computation 
of the Fokker-Planck equation, as a method for stellar dynamics simulations. Additional 
benefits of this approach are that the Fokker-Planck equation takes a particularly clean 
form and that the procedure of orbit-averaging is then a conceptually simple averaging of 
quantities over the 2tt change in (an) orbital angle variable, as will be seen in Chapter [3l 

With this form it is also possible to straightforwardly incorporate terms describing collisional 
stellar mergers into the Fokker-Planck equation. The main drawback is that a large number 
of conversions are required between the {E, J, J^) coordinates in which dynamical quantities 
are calculated and the (I, J) space in which the system evolution is tracked by the Fokker- 
Planck equation for /(I, J). Integration of the two-dimensional Fokker-Planck equation - 
even in its more clean form, as when expressed in term of actions / and J - is sufficiently 
challenging that this tradeoff remains beneficial. Effectively the additional complexity of 
the dynamics has been separated out of the Fokker-Planck equation itself and into these 
coordinate conversions. The result is a three-step simulation process. First the Fokker- 
Planck equation is used to update the distribution function /(I, J) at a new timestep. Then 
the effects of any additional physical processes - in this case, collisional stellar mergers - are 
calculated and /(I, J) updated accordingly. Finally, given the new /(I, J), consistent values 
for physical quantities such as the stellar density and gravitational potential are found in 
terms of either E and J or the spatial coordinates. 

1.3 Overview of Resulting Model and of Remaining Chapters 

Given all the above considerations, the basic form of the model developed for this thesis 
has the following characeristics: 

1. The system to be studied is a fairly dense stellar cluster, using a generalization of the 
initial conditions from Quinlan and Shapiro [7], with an observationally-motivated 
initial mass function for the stars which make up the cluster; 

2. Unlike most prior work on dense clusters, an overall rotation of the system is included, 
with a value determined by observation; 

3. In order to allow for the inclusion of desired additional physical effects beyond the 
general gravitional interactions of individual stars in the cluster - i.e., a stellar bar 
perturbation, and collisional merging of individual stars - the Fokker-Planck equation 
is employed to track the system's evolution; and 
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4. To capture the physics required while still keeping the problem numerically tractable 
the Fokker-Planck model is orbit-averaged, and any third integral neglected, so that 
the system is represented in two dimensions. The canonical actions (/, J) are employed 
as integrals of the motion in order for the above-mentioned extra physical effects to 
be incorporated straightforwardly. 

With this setup for the model, the basic question addressed here is under what circumstances 
- using which physical effects, starting with what initial conditions - can the system produce 
a star which can serve as seed object for an intermediate-mass black hole and eventually a 
supermassive black hole. For this to be possible, the star must have sufficient mass, i.e., 
25OM0 or more, and it must be formed before the end of the main-sequence lifetime of 
large stars, i.e., within the first 3 Myr or so of the system's evolution. The simulations will 
show that a bar-like perturbation can effect sufficient angular-momentum transfer to allow 
the creation of a likely seed object in most rotating clusters whose density and velocity 
dispersion match what is seen in the centers of galaxies, but only if the initial mass function 
is also set to match that observed in dense clusters - which is more top-heavy that the 
general field-star IMF. 

The remaining chapters are organized as follows: 

• Chapter 2 outlines the dynamical and gravitational-potential physics required for the 
Fokker-Planck model, lists the numerical considerations needed to compute it, and 
outlines the set of possible initial conditions. 

• Chapter 3 derives the diffusion coefficients that incorporate the gravitational effects 
of the system's individual stars and of any stellar bar perturbation into the Fokker- 
Planck equation for the evolution of the system, as well as the population effects 
of collisional mergers. The differencing scheme employed to compute the model is 
developed, and a brief consideration of the heating and merging effects of binary stars 
is given. 

• Chapter 4 describes the tests performed on the various aspects of the model, both to 
verify its numerical and physical validity and to establish reasonable values for any 
numerical input parameters required. 

• Chapter 5 details the results of the simulations for the various possible choices of 
initial conditions and of which physical effects to include. 

• Chapter 6 then discusses the simulations' results from an astrophysical standpoint, 
gives comparisons to what would be expected from previous studies and from general 
considerations such as timescale arguments, and looks forward to what upcoming 
observations and future models may be able to provide. 

• Tables of frequently-used symbols and subscripts are given in the Appendix. 

In the end, the simulations show that dense rotating stellar clusters are a feasible path to 
the production of seed objects for massive black holes in galaxies. All simulations that were 

performed with a top-heavy IMF - such as is observed in dense stellar clusters - resulted in 
at least one sufficiently massive star (and usually several) being formed. The result holds 
for a range of initial cluster configurations which were chosen to match typical dense nuclear 
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clusters and galaxy centers. Only for models representing giant elliptical galaxy cores was 
a massive object not seen before the largest stars evolved off the main sequence. For these 
systems one would need to rely either on post-main-sequence dynamics within the core, or 
on the presence of a dense nuclear stellar cluster which could evolve independently. 
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Chapter 2 



Dynamical 
Aspects of 



and Gravitational 
the Model 



2.1 Overview 

The first two sections of tliis chapter describe the simulation model's dynamical framework, 
specifically the calculation of orbital parameters e.g., endpoints, frequencies etc. These will 
be needed in Chapter [3] to determine the effect of dynamical friction on the ensemble of 
stars in the cluster, as described respectively by the diffusion coefficients and the stellar 
distribution function in the Fokker-Planck equation. The remainder of the chapter then 
develops the calculations of the cluster's physical properties, such as the mass density profile 
p{r) and the gravitational potential $(r), and ends with a summary of the possible initial 
conditions used for these quantities. 

The fundamental units employed in the model are km/s, Mq, and pc. With these choices the 
unit of time is lpc/(km/s) = 0.9778 x lO^yr ~ IMyr, and Newton's constant of gravitation 
is G= 4.299 X lO'^ (km/s)^ pc/M©. 

2.2 Dynamics 

The most convenient way to parameterize stellar orbits is to use to the canonical radial and 
angular actions, which together will be referred to as the action vector Ij, with components 
Ii (the radial action) and I2 (the angular action). I2 is equal to the traditional angular 
momentum J, while Ii can be defined from the radial momentum: 



in which orbital endpoints rp, are the star's periapse and apoapse distances from the 
cluster center, E is the star's total orbital energy and $(r, 0) is the gravitational potential 
at distance r and polar angle 9 from the center!^ Azimuthal symmetry of the cluster's main 

^Shlosman |28) points out that the Jacobi energy E of (|2.ip isn't a true integral of motion, even for 
individual stars, when there are nested bar perturbations: the underlying potential is then not a constant. 




(2.1) 
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stellar distribution (and so also of the potential) is assumed. 

The advantage of using the canonical coordinates is that the stellar distribution function, 
when expressed in terms of these actions, is invariant to adiabatic changes in the potential; 
hence we do not have to solve for a new distribution function just because the potential 
evolves over time |42j . As seen above however it is often necessary to deal with the tradi- 
tional variables r, 9, (j), E, J, when calculating orbital parameters or any physical quantity 
that depends on the value of the potential. Alternatively, when the actions (Ii,/2) are 
known, (12. ip implicitly defines the stellar orbital energy E. 

Following Tremaine and Weinberg [30], we can move between action space and position 
space via the relations that follow. The radial and angular orbital frequencies are found 
similarly to Ii : 



1 1 dr 



ni TT J^^ [2(S- «>) - j2/r2]V2 
^2 J dr 



-K Jr^ r2[2(^-$)- j2/r2]l/2 

From these we can relate the dynamical quantities in the two coordinate systems: 



( dl\ _ J_ fdl\ ^2 



\dE J J Qi' \dJ p 



(2.2) 
(2.3) 

(2.4) 



dE 



Qj (2.5) 



where the last relation generalizes to the case including a third integral of motion if we 
define I3 = Jz and note that ^3 = 0. 

Although orbit-averaging will remove any direct dependence of the diffusion coefficients 
on the canonical angles wi and W2, the angles are needed for intermediate steps in the 
coefficients' calculation. (The angles are also needed for some of the tests of the coefficients' 
validity.) They are given by 



c 



and 



{Sl2-J/rHdr\ , , 



C 



The azimuthal angle is not of interest except in that it contains a random initial phase 

As shown in Fig. 12.11 is the angle swept out by the star's position vector as measured 
from the point where it crossed the 6 = tt/2 plane: 

sin-f/'sin/? = cos6', (2.8) 



This is expected to be the case for clusters such as those modeled here when near the center of a galaxy 
[36| : however, the inner bars decouple dynamically from the larger-scale galactic bars [37] which will vary 
on much longer timescales, and so this nicety is expected to produce only a secular change in the underlying 
potential that averages out over long timescales. 
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Image source: Tremaine and Weinberg jSOj . 
Figure 2.1: Euler angles for the orbit of a given perturbing potential. 

in which (3 is the orbital inclination, i.e., Jz = Jcos/3. Integration path C is along the 
orbit, i.e., from to and back to when traversing an entire orbit. For our purposes 
the quantity (^2 — 4') is what will be needed most often. The above relations are intuitive 
given that [2{E — $) — J'^/r^]^/^ is the radial component of the velocity and the canonical 
angles obey (fwj/dt^ = 0. 

2.2.1 Numerically Solving for Orbital Endpoints 

The general method for finding the orbital endpoints r^, is to solve numerically for the 
roots of = ^{E — <!>) — J^/r^] = 0. For small values of either action component Ii or J 
this may not be numerically robust and other techniques must be used. If {E — $)/|$| is 
miniscule {i.e., < 10^^), the star is considered to be "fixed" at the origin. For (i? — $)/|$| < 
10~^ a quadratic approximation to the potential is employed and the endpoints solved for 
analytically. The endpoint-solver takes as its natural inputs E and J, but it also requires at 
least a reasonable guess of Ii so it can distinguish between more-radial and more-circular 
orbits. 

For very-nearly-radial orbits (taken to be |J| < ^^toi^)) '''p is set to a value of 10~^ in 
order to avoid numerical problems at exactly r = 0. □ The input parameter Ztoi = 10~^ 

^Alternatively, for nearly-radial orbits Ii/\J\ < 10"'^, an epicyclic approximation for finding rp and Va as 
well as Q.j was tried. This was not found to be an improvement and results with the epicycle approximation 
are not reported here. However, for numerically tiny values of the components of the action (typically 
-^j ^ 5 X 10"'*) I.e., for stellar orbits restricted to the cluster center, an epicylcic approximation 0^ is 
employed: {yi2f — d^$/dr^\r=o and £7i = 2Q2, although in practice no values of 7, on the grids used trigger 
this case. 
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typically, although other values were tested to make sure the results were robust against 
such a change. 

For moderately-circular orbits (those with Ii/\J\ < 0.17), first a value of r for which Vr > 
is found; then and are determined by searching out in both directions until Vr changes 
sign. For more radial orbits, the entire allowed range of r is searched for subranges over 
which Vr changes sign. In both cases, the bracketing of the ranges over which Vr changes 
sign is then made increasingly finer until convergence is reached. 

2.2.2 Calculating Orbital Frequencies and other Dynamic Quantities 

Once orbital endpoints and are known then ()2.3p and (j2.2p can be used to find orbital 
frequencies directly (and similarly (12. ip gives the action components Ij if they are not 
yet determined). In practice doing integrals is very expensive in processor time; thus in 
general ilj is calculated from (12. 5p . When values for dE/dlj are not yet known for a given 
simulation timestep, however, ()2.3p and ()2.2p are employed. 

Just as in the orbital endpoint calculation, boundary cases require special treatment. When 
the orbit is highly radial (taken to be when J/Ii ^ 5 x 10"'*), ()2.3p becomes unstable and 
so the limiting case of ^2/^ = ^ is used. Finally, for orbits that are very nearly circular 
(Ii/J < 5 X 10~^), even if the solution for rp and is obtainable, doing further calculus 
such as (j2.3p and (j2.2p is numerically difficult and so in this case the orbit is taken to 
be exactly circular with Tp = = r^, the radius at which dv^/dr = 0. The above Ii/J 
requirement is sufficiently strict that this will not affect any dependent calculations. 



2.3 Calculation of the Density 

2.3.1 Using the Distribution Function in (E, J^) Space 

Traditionally, the (smoothed) stellar mass density p is given by the "full" 6-dimensional 
integral of distribution function F over all velocity space - which gives the total number 
density of stars - multiplied by the per-star mass m: 

p{r) = m j d^vF{r,v) (2.9) 

which results in a total mass M for the cluster of: 

M = j p{r)dV = m j d^r d^v F{r,v). (2.10) 

Previous studies assumed a stellar cluster with spherical symmetry and used {E, J^) as the 
coordinates of choice; using the transformation of Cohn [44j, in {E,J'^) space the above 
expression becomes 

M = 8TT^m j d(j2) j dEQ.:y^F{E,J). (2.11) 

With azimuthal symmetry, 27r of the above comes from integrating over the velocity-space 
azimuthal angle 4>y. Expanding d{j'^) and applying (for fixed J) dE = Q,idl from (12. 3p . the 
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total cluster mass when integrated over the actions is thus 

M = 16-71^771 j JdJ j dl F{E{I,J),J). (2.12) 

2.3.2 Using the (E, J^)-space Distribution Function in Action Space 

In order to take advantage of the invariance of the actions, we define a distribution function 
analogous to F but expressed in action space: 

M = m j dIdJf{I,J). (2.13) 

Note that this is by construction an "orbit-averaged" distribution function ([42j, [44j ) and 
so all angle-dependence has already been averaged out. Comparing with ()2.12p one obtains 

/ = IQtt^JF (2.14) 

as the conversion between the action-space / and the earlier F. 

To calculate the spatial density p, one can either perform the integral of (j2.9p directly, or 
first convert it to a 2-dimensional orbit-averaged form that uses /, similar to (j2.13p . The 
former requires use of the velocity volume element (see e.g., Cohn and Kulsrud [47j) 

2^ ,3 , JdJdE 

d^v = Att — K (2.15) 

and so, if using the two-dimensional / of (j2.13p and assuming spherical symmetry, 

p{r) = 47rm / dE / 1^F{E,J) (2.16) 
J$ Jo r^Vr 



Anm / JdJ / dI—^F{EiI,J),J) 
Jo J r Vr 



where the endpoints of the / integral are those energetically allowed. This has the advantage 
of being a direct method of calculation, but is otherwise undesirable for this study for several 
reasons: the l/r"^ dependence means that the p calculation is less robust close to the cluster 
center; the singularities at the orbit endpoints due to the 1/vr dependence are particularly 
difficult to deal with in the dIdJ form of (j2.16p , and spherical symmetry is not an intended 
assumption for this model. 

2.3.3 Calculating the Density in Action Space 

In order to deal with the above issues it is necessary to develop a differential in terms of 
the action coordinates themselves. Expanding d^v in spherical coordinates we have 

d'^v = v'^dv d cos By d(f>y. (2-17) 
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Defining the "tangential" velocity Vq+v'^ = v'^ = {J/rY, then dJ = rdvx = rvdcosO^ and 
so 

d^v = 2vdvd(t)v— (2.18) 
r 

and 

d^v = ^'Kvdv — . (2.19) 

Note that 2tt in the above again comes from integrating over azimuthal angle and an 
extra factor of 2 has been inserted in order to allow use of the usual convention that J > 
only. Finally, with Vtidl = dE = vdv, 

dJ 

d^v = ATT^idI—. (2.20) 

'i,=0 ^ 

Inserting this in ()2.9p and using the conversion from F to f, 

m f didj ^ „, ^ , ^ 

where as always, the action integral endpoints are determined by which orbits are energeti- 
cally allowed at r; this is indicated by subscript "r" on the integral sign. This form for the 
calculation of p is more suited for the overall action-space-based model of this study than 
is (j2.16p . and is what will be used throughout; the cost of this choice is that the orbital 
endpoints are now more complicated to calculate. 

In what follows, it will often be convenient to use the notation (/i,/2) in place of (/, J). 
Generalizing the expression for p to the case in which there are multiple values for the stellar 
mass, 

Pir) = ^Jf-^ n,Y.m,U{h,h). (2.22) 

where subscript q refers to a particular stellar species {e.g., if the lowest-mass stars in the 
simulation are of solar mass, then mi = IMq and those stars' distribution function is /i). 
The ranges of integration of both of the Ij are the ranges of energetically allowed values, 
and so depend on E{Ii,l2) and on <I>(r). 



2.4 Rotational Velocity and Orbital Inclination 

It is possible to deduce an expression for the rotational velocity zu^r) of a given stellar orbit 
using only the velocity dispersion and the assumption of velocity isotropy. However, such 
a relation would depend more heavily on the assumption of isotropy than does anything 
else in this study, and in the end it does not yield any insight into the distribution of orbits 
over inclination angle /3 (where cos/3 = J^/J), which will also be needed later. Instead, we 
implement a procedure for first finding the distribution over /3 and use that to determine 
zu. 
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2.4.1 Orbital Inclination 



Even in a spherically-symmetric gravitational potential, with an overall rotational compo- 
nent to the stellar velocity distribution some quantities will depend on the resulting distri- 
bution of orbital inclination. It is convenient to parameterize the inclination dependence of 
/ as having factors that are even and odd in a = ^ — /3: 

h{a) = g{\a\) + e{a) (2.23) 

where j^^^i^dah{a) = jZ^^i2dci g{\a\) = 1 and 0(— a) = —©(a). Thus defined, h{a{j3)) 
can then be applied as a weighting function in averaging any /3-dependent quantity, which 
is essentially a way of defining an effective distribution function: 

/^(/, J; P) = h{a)f{I, J) = [g{\a\) + Q{a)]f{I, J). (2.24) 



Without any prior knowledge of the rotational state of a stellar system other than the net 
amount of rotation, DeJonghe (|27j. |48j ) has shown that the statistically most probable 
form for Q{a) is 

e(a) = tanh f^JM] g(\a\) (2.25) 



where 6 is a parameter determined by the total rotational angular momentum. Other than 
being an even function of its argument, so far 5(|a|) is unconstrained. A particularly useful 
choice is to assume that for a given J, all possible values < | J^l < J are equally likely |30j . 
consistent with the assumption of isotropy. This is equivalent to setting ^daj) = icosa : 



Ha) = - 



1 + tanh 



bJ sin a 



cos Q 



(2.26) 



or, rewritten in terms of inclination angle /3, 



1 + tanh 



bJ cos /3 



sin /3. 



(2.27) 



It remains to determine the value of b. For notational convenience we define mass-weighted 
distribution functions, i.e., f = Yliq'^^qfq similarly F = Ylq'^q^q- Then we define the 
"population" of a given orbital energy as 

P{E) = / dJF{E,J) (2.28) 

'-^ min 

and the total angular momentum for a given energy as 



'max 



Jtot(.E)= / dJJF{E,J). (2.29) 

-^min 

For a distribution that has no J-dependence (and hence, no rotational streaming), the 
expected total angular momentum is simply 

1{E) = ^ (J^^ + Jmin)P{E) (2.30) 
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and so one can deduce the excess of angular momentum at a given E: 



AJiE) = JtotiE)-J{E) 



(2.31) 



so that AJ{E) = for no net rotation. This allows for the following scheme for determining 
b: for a given orbital energy E, choose b{E) such that the resulting amount of net rotation 
is consistent with AJ{E): 



AJ{E) 



dJF{E,J) / daJ^h{a;E) 

TT 

2 / dJJF{E,J) / ' dasmaQ{a;E) 



dJJF{E,J) / (i(sin a) sin Q tanh 



bJ sin a 



dJJF{E,J)[—j 



2 \ 2 i-bJ/2 



dxx tanhx 



which can be expanded to read 



dJF{E, J) 



J 



2 \ 2 i-bJ/2 

bJ 



dx X tanh x 



0. 



(2.32) 



(2.33) 



Root-finding over the above integral equation is a computationally expensive process; before 
resorting to that, we first attempt an iterative solution by expanding xtanhx for small 
a; < I and solving the resultant equation for the leading term in 6, using the previous 
iteration's b value in the higher-order terms. If this does not converge we do rootfind on 
the exact expression above. For larger values of &J/2, i.e., bJ/2 > ^ the integral over x 
itself asymptotically approaches ^, simplifying the process. In all cases the resulting b{E) 
is then verified to satisfy (j2.32p : with this 6, (j2.27p then describes the orbital inclination 
distribution of orbits with energy E. 



2.4.2 The Mean Rotational Velocity 



Further defining the population p at a given r as 

p{r) = f dIdJf{I,J) (2.34) 

Jr 

finding the mean rotational velocity w at a given r is then straightforward: 

1 f f'^ J 

tu(r) = —- dIdJf{I,J)lda — h{a{E{I,J))) (2.35) 
p[r) 7_| r 

1 r ~ / 2 \ ^ i'^-^/"^ 

= — — dIdJJf(I,J){ — ] / dxx tanh X (2.36) 
rp{r)Jr \bjj Jo 

in which b = b{E{I, J), J) as found in ^2.4.11 Descriptively, w{r) is the mean speed at which 
stars stream in the azimuthal direction at a given distance from the system center, and is 
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found by averaging J^/r over all inclination angles /3 as weighted by the relative density 
2.4.3 The Coordinate Grid 

The self-consistent gravitational potential <I>(r) needs to be updated after each time step 
in the simulation, through the iterative solving of Poisson's equation. Up to this point, 
the dynamics - including the calculation of the orbital inclination distribution h{/3) - have 
assumed a spherical gravitational potential <I>(r) and (smoothed) density p{r), i.e., that the 
cluster's overall ellipticity e = 0. This is done simply to allow the construction of a self- 
consistent model using the two-dimensional stellar distribution function /(/, J), as adding 
a third integral would be computationally prohibitive. 

With the present definitions of dynamical quantities and construction of /i(/3), setting e = 
is mathematically proper [27] and yields a self-consistent model. It is possible that e 7^ 
may be more physically motivated, although N-body studies do show that "spherical stellar 
systems can rotate very rapidly without becoming oblate" [23j. In order to see how much 
effect the assumption (by construction) of sphericity has, we can add in allowance for 
ellipsoidal isodensity surfaces in an ad hoc manner, using the presumption that the first- 
order effects of e(r) > on bulk properties {e.g., anything that depends directly on p or 
on <I>) will outweigh the (small, for small e) errors it introduces into the dynamics properU 
Thus unless otherwise specified, ellipticity is pre-defined as e = for all results - although 
provision for the future consideration of the case with a preset e(r) 7^ 0, which may be 
allowed to vary along with the coordinate grid values, is built into the model and is the 
basis for the development presented in ^2.51 of the gravitational potential and in Appendix 
IA.2l of numerical aspects of the model. 

Thus in the general case there will be a ellipticiy e(a) associated with a given isodensity 
surface at ellipsoidal radius a; this must be taken into account when calculating quantities 
that depend on the local position]^ Equating the exterior value of the potential (i<I>e(ffl) 
of an ellipsoidal shell of mass dM with the analogous value d^Qir) for a spherical shell of 
identical mass one obtains 

GdMsin-^e_ G dM 
ae r 

( sin~^ e\ ,^ 1 On 
a = r ~ r(l H — e ) 

V e ; 6 

where the approximation is used throughout for very small (< 2 x 10^'^) values of e. From 
here on, r can be taken as a "dummy" variable, used for convenience to label gridpoints; the 
shell's semimajor axis length a from (j2.38p is the physically-meaningful radial coordinate. 

In cylindrical coordinates, c? = V? ^ \+e{a) ' '^^^ ^ ^I'id is allowed to dynamically update 
with each time step, meaning that the values of a} chosen to be grid points can change at 
each timestep so that regions of larger d^jdc? and dpjdc? are given a greater density of 
grid points [39]. When e = then a = r; this is the case for the majority of simulations 

^The only inherently non-spherical aspect of the model's potential or density, namely a "bar perturbation" 
in the stellar distribution's potential, will be considered in ii3.6.3l 

*Note that, despite allowing for a net overall rotation in the stellar population, the procedure for calcu- 
lating the non-constant distribution of orbital inclinations h{P) does not alter the overall energy balance. 
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(2.37) 
(2.38) 



studied here, and for those the coordinate grid is referred to as the "r^ grid" and not as the 
"a^ grid" . The numerical consequences of having a finite number of coordinate gridpoints 
are discussed in ^4.7.2[ 

2.4.4 The Velocity Dispersion 

The one-dimensional root-mean-squared ( "rms") stellar velocity dispersion cjo at a given 
radius a or r is found similarly to w: 

al = (1 - 5)al = ^1-^ j^dldjj{l, JME - ($) J - (4>J (2.39) 

where 5 is an anisotropy parameter that allows for a difference between the radial velocity 
dispersion ui, and (Tq- Subscript "a" on bracket {■) ^ denotes the possible averaging of the 
contained quantity over the surface of the isodensity ellipsoid with semimajor axis a. To 
wit, note that ()2.39p assumes when 5 = 0, do is identical to ai = {v^.)^/"^. This is necessary 
given that the restriction of working with a two-dimensional /(/i,/2) means there is not 
full information on the three-dimensional velocity distribution, and so we must rely on an 
assumption regarding the degree of velocity anisotropy in order to calculate a^. It is also 
why it is not advisable to depend on this calculation of al when determining either the 
distribution of inclinations h{(3) or the amount of stellar rotation; those quantities are of 
central importance to the overall model whereas do is less so. 

The expression ()2.39p gives the velocity dispersion for all stellar mass species considered 
together, but it is just as easy to perform the same calculation for a particular stellar mass 
value by restricting the sums in p{r) and / to that single species. 

2.5 Gravitational Potential 

As the simulation progresses, new values for the stellar density distribution and the grav- 
itational potential must be determined at each timestep in turn. The new potential at 
each point on the radial-coordinate grid at a given timestep can be calculated in a manner 
analogous to that used by Cohn ^44j- We write the (inverted) Poisson equation 

^new ^ £"lp[ynew. ^new] (2.40) 

where is the inverted Laplacian, derived below. The "new" superscript denotes that, 
starting with a newly-updated stellar distribution function /, (|2.22p and then (j2.40p are 
used to find self-consistent values for p and Note that the allowed ranges of Ii and I2 in 
(j2.22p depend on ^{R,z), hence the implicit dependence of p on <^ in (j2.40p . For the first 
step in the interation = is typically used in the right hand side of (I2.40p : if this 
fails to converge a guess for ^^"^^ based on a comparison of previous timesteps' potential 
functions is also attempted. 

As done by Cohn [2] 5 we employ the "Aitken (5^" process [50] to accelerate convergence. 
This procedure is iterated until convergence is achieved, typically to an accuracy of ~ 2% 

'^For {v^)^, simple geometry yields = 3^ + 3^:^(1,^:^) = + 3(i3^]- The average value of the 

gravitational potential {$)^ over the ellipsoid surface must be calculated numerically. 
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for every point in the Aitken-processed grid. A final check of conservation of overah 
energy is also made. 

When the gravitational potential $ is taken to be spherically symmetric - as is the case 
for most of the simulations performed - the procedure is conceptually straightforward, if 
numerically complex. The mass density is updated from the new distribution 

function using (j2.22p and then the solution of the Poission equation is simply 



after which <I>'^''™] and <|)°™ are iteratively solved for until convergence is achieved. 

For the more general, but rarely needed case in which there is an overall and variable 
ellipticity e(a) to the potential, a somewhat more complex procedure than given above is 
required. One still starts with the inverted Poisson equation as given in (j2.40p but now 
the derivation of the inverted Laplacian is more involved. Details are given in the 
Appendix. 

2.6 Initial Conditions 

It is standard in stellar-dynamical work to first consider a Plummer sphere distribution 
|44j before moving on to other possible potential-density pairs. The Plummer sphere is 
the reference potential for this study: it fits the light curves of globular clusters [5Tj. An 
alternative choice of initial potential which fits the surface brightness of galactic spheroids 
is the "7 = 0" model [52j. These two choices have the shared advantage of providing well- 
defined values for the initial distribution function /(/, J) in addition to the potential $(r) 
and density 

To introduce rotation into either of these density-potential pairs requires an ad hoc al- 
teration of the distribution function as described below in ^2.6.21 which is followed by 
descriptions of the stellar bar model and the possible choices for the stellar initial mass 
fuction (IMF). 

2.6.1 Potential-Density Pairs 
Base Model: The Plummer Sphere 

Previous studies of dense-cluster dynamics have often started with a Plummer sphere ([7], 
|23j . |53j . |41j). going as far as refering to nonrotating single- mass Plummer spheres as a 
"standard testbed" [Mj- Following Quinlan and Shapiro [T], the gravitational potential and 
stellar mass density for a Plummer model of total mass M are 




(2.41) 



$(r) 



GM 



(2.42) 



V ~r ' core/ 



and 



3M 



(l+rVri,J-| 



(2.43) 



core 
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The corresponding energy-space distribution function for mass species niq is 



F(E)- ^^^^g ^7/2 (2 u) 

where Nq denotes the total number of stars of individual mass niq. Note that a general 
non-rotating distribution function is not necessarily independent of J, but the Plummer 
model's is. 

In particular the specific models first employed by Quinlan and Shapiro [1], such as those 
listed in Table \27l\ have also been commonly used {e.g., as done by Rasio et al. p^) for 
later studies of nonrotating systems, and so are modified for use here as described below in 
§2.6.21 The models in Table 12.11 have two desirable properties: the initial relaxation time 
tr < 1 X 10'^ Myr so the total collapse time is expected to be < 1 x 10"^ Myr; and the 
escape velocity satisfies o"i < fesc — 600km/s and so stellar collisions result in coalescence 
of the bulk of the combined stellar material, not disruption and dissipation of it [7|. The 
naming convention is that models with the same middle digit in their name share a common 
value for initial velocity dispersion cro(O), while the final letter indicates the initial central 
relaxation time tr{0). 

Relating the models to astronomical systems, model ElB best describes globular clusters 
and is not studied extensively here. E2A and E2B could be very-dense globular clusters 
[55], or the nuclei of bulgeless spiral and dwarf elliptical galaxies (as shown in Preitag 
et al. [56], [57]). The densities of models E4A and E4B match those of massive young 
stellar clusters [56] • The cores of giant elliptical galaxies are observed to have line-of-sight 
velocity dispersions o"los near 350 km/s [58j; this is not quite a match for the do = 400 
km/s of models E4A and E4B. However, near the observed centers of galaxies (Tlos is more 
directly comparable to o"i (as given bv I2.39P and it will be shown in ^5.2.21 that for model 
E4B, (Ti(0) ~ 340 km/s, a good match with the observations of giant elliptical galaxy cores. 
Previous simulations of nonrotating systems have found that a very massive star - identified 
as a possible precursor to an intermediate-mass black hole - can form in a system with a 
velocity dispersion of "many hundreds of km/s" [19], so the high value for the velocity 
dispersion should not be an impediment. 



Alternative Model: the "7 = 0" Sphere 

The potential-density pair for the "7 = 0" sphere is 

(2.45) 



^>^(r) 



2r, 



(r + r, 



,2 

core / 



p^{r) = ^^£2£^ (2.46) 



with distribution function 



3MA^„ /, ,i3-4e , / 2e \ Ercore , 
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Of other possible models, note that the potential of Jaffe's [59] model, and the density of 
Hernquist's [60] model, diverge as r — )• and so are inappropriate - or at least inconvenient 
- for use here. 

Carrying over values of M and Tcore from Table 12.11 to the case of the 7 = sphere, only 
one resulting model has central density p{Q) and velocity dispersion ai in the required 
ranges: "G2A", i.e., a 7 = sphere with the same values for Af and rcore as model E2A. If 
instead of equating rcore one uses common values of the half-mass radius ri/21 a potential 
model "G3C" results with core radius rcore = 0.20 pc. Model G3C has the same total 
mass M = 1.1 X IQ^Mq and the same half-mass radius ri/2 as model E2B, albeit with a 
very large initial central density of /9(0) = 1.4 x IO^Mq/pc^; this model is likely of most 
use as a demonstration case rather than as a representation of a realistic astronomical 
system. However an intermediate model with rcore = 

0.29pc {i.e., half that of model E2B's 
?"core = 0.58pc) features a less extreme initial density p(0) = 5.0 x IO^Mq/pc^. This model 
has an initial relaxation time similar to that of model G2A but a somewhat larger initial 
central velocity dispersion and is thus labelled "G3A". 

Numerically calculated values for velocity dispersion cri, density p and relaxation time ty. 
for these 7 = sphere models are given with the overall simulation results in §5.31 



2.6.2 Introducing Rotation 

As described in §1.1. 3^ rotation due to cosmological tidal torques is ubiquitous in the uni- 
verse and is only weakly dependent on how a particular system was formed. Barnes and 
Efstathiou [2T] considered various models for formation of objects in the early universe, and 
in terms of Peebles' dimensionless spin parameter 

\ _ -^rot I Wgrav | ^ .„s 

they determined that the typical value due to cosmological tidal torques is A ~ 0.05, almost 
independent of the perturbation spectrum. At the stellar cluster scale required here this is 
consistent with the measurement by Cervantes-Sodi et al. |61j . who found for a sample of 
SDSS galaxies that A = .04 it .005 with a weak trend of increasing A with decreasing mass. 

In order to introduce rotation into the simulations, we alter the distribution of Fq{E) so 
that the total angular momentum of the cluster Jrot produces the desired value for A, while 
conserving Nq. This results in a "tilted" F{E, J) in which the resulting excess total amount 
of J is attributed to Jrot^ 

with the tilted distribution function 

Fq{E,J) = {l + Q,ot)Fq{E) (2.50) 
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Model 


P(0) 


fTo(O) 


M 


'"core 


t,(0) 


Corresponding 




[Mq/pc^] 


[km/s] 


[Mq] 


[pc] 


[yA 


Astronomical object 


E4A 


3.0 X 10« 


400 


1.8 X 10^ 


0.24 


4.6 X 10^ 


core of a 


E4B 


1.0 X 10^ 


400 


3.1 X lO'^ 


0.42 


1.4 X 10^ 


giant elliptical galaxy 


E2A 


4.0 X 10'' 


200 


6.2 X 10*^ 


0.33 


4.6 X 10'' 


nucleus of dwarf elliptical 


E2B 


1.3 X 10^ 


200 


1.1 X 10^ 


0.58 


1.4 X 10^ 


or bulgeless spiral galaxy 


ElB 


1.8 X 10*^ 


100 


3.6 X 10*^ 


0.79 


1.4 X 10^ 


globular cluster 



Table 2.1: The set of initial condition models. For purposes of the numerical simulation, 
the quantities taken as fundamental parameters are total mass M and core radius rcore- 
The central density /o(0), velocity dispersion cTo(O) and relaxation time tr-(O) listed here are 
the resulting analytic values for a Plummer sphere distribution. Model names are taken 
from Quinlan and Shapiro [Ij. 

in which the amount of tilt QrotiE, J) is given by 

In short, this procedure simply replaces Fq{E) with a new Fq{E, J) that conserves number 
and gives the same total angular momentum, but that has a non-flat J-dependence of the 
form (J — (J)^)^ which results in a net total rotational angular momentum Jrot- 

Because the (/, J) grid (on which the Fokker-Planck coefficients are calculated) is square 
but the corresponding (E, J) grid (representing the space on which dynamical quantities 
must be calculated) is not, we must take care in how A J is determined in the rotating case. 
In the above, {J) ^ = ^(Jhi(-E') + Jio{E)) is the average value of J on the grid for a given 
E. Similarly, A J is the span of J values over which Fq(E) is being tilted; to strictly satisfy 
(I2.49P one would take it to equal Jhi{E) + Jio{E) but in practice this results in values of 
E which do not span a large range of J being "overloaded" with more than their share of 
rotation. Thus insteacifl we simply set A J = (Jmax — Jmm)- Polynomial power C > 1 is a 
free parameter that sets the shape of the tilt. We usually apply a linear tilt C = 1- 

The new rotating distribution F{E, J) is no longer a solution of the Poisson equation that 
matches the Plummer potential and density profiles given by (|2.42p and (|2.43p . This "prob- 
lem" is easily overcome by using the mechanisms of the main simulation itself, as described 
earlier in this chapter, to find the <I>(r) and p{r) that do correspond to the new F{E, J), 
before starting the simulation proper based on these initial conditions. 

2.6.3 Initial Mass Function 

The spectrum of masses which is input to the model is expected to play a large role in 
its dynamical evolution [54] . A simple power-law dN/dm oc m^" represents perhaps the 

^The cost of this change is that the resulting value for A of the newly-determined distribution F{E, J) 
is slightly lower than intended. This would happen in any case however, as the new distribution also has a 
lower value for |Wgrav| than the original Fq{E) did, due to the additional overall rotation. All values for A 
quoted in Chapter [5] are "true" values in that these two effects have been taken into account by recalculating 
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simplest form for the initial mass function (IMF). The traditional Salpeter IMF has a = 
2.35. Other observationally-determined IMFs include the Miller-Scalo |5l] : 



dN 
dm 



m"^-^, 0.1 < m/Mo < 1.0 

1.0 < m/Mo < 10 (2.52) 
10 < m/MQ 



and the Kroupa IMF [62]: 




0.08 < m/MQ < 0.5 

0.5 < m/MQ < 1.0 (2.53) 
1.0 < m/MQ 

A problem with the above IMFs is that none was determined for the specific case of dense 
clusters near centers of galaxies. However conditions in the centers of galaxies are sufficiently 
different from those elsewhere that all but the highest-density molecular clouds (number 
density n > lO'^/cm^) will be shredded by tidal forces p^. Recent observational work on 
such IMFs - most notably for the Arches cluster located 30pc from the center of the Milky 
Way, as well as the Quintuplet and Central clusters, and R136 in the Large Magellanic 
Cloud, has led to the following picture of the IMF in the centers of galaxies (Figer |12) . 
[63], Oey and Clarke [64J, and references therein, but see Elmegreen [65] for a dissenting 
opinion) : 

• galaxy-center IMFs are top-heavy, with a ~ 1.7 — 1.8 for M > (6 — 10)Mq, or ~ 
1.8 — 1.9 when differential extinction is taken into account; 

• stellar masses below ~ 2Mq do not add appreciably to the total mass of the stellar 
cluster; 

• there is a fundamental, initial upper mass cutoff of ~ I5OM0; 



central densities can reach ~ IO^Mq/pc^ 



Notable exceptions to the high-mass cutoff include the Pistol Star and its "twin" FMM362 
in the Quintuplet cluster; despite a cluster age of at least 4 Myr these stars have masses in 
the 150 — 2OOM0 range with an expected stellar lifetime of 2.5 — 3 Myr. One explanation 
is that these stars are the product of stellar mergers of stars in the IOOMq range |12j . 



Thus an IMF for the Arches cluster, representative of dense clusters near the centers of 
galaxies, can be expressed as 

dN j m-\ 2<m/MQ<8±2 

^ « I ^-i.8±o.i^ 8 i 2 < m/MQ < 150 ^^■^^> 
and a simplified version using somewhat overly-conservative choices of parameters is 

^ oc 2 < m/MQ < 150. (2.55) 

dm 

It is the above IMF that is used for the "Arches" simulations in this study. Some other 
cases here use the Salpeter IMF of a = 2.35 but most start with the Kroupa IMF which for 
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a lower-mass bound of IMq corresponds simply to a = 2.7; both of these were also used 
by Amaro-Seoane [H] for an analogous study of dense but nonrotating clusters which had 
already formed massive central objects, i.e., black holes or supermassive stars. 

With all relevent dynamical quantities now calculable, the Fokker-Planck equation's coeffi- 
cients for the system can now be determined; the validity of the overall simulation technique 
and its various components tested; and the actual simulations performed using the initial 
conditions given above. These steps comprise the next three chapters respectively. 
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Chapter 3 



Derivation of the Fokker-Planck 
Diffusion Coefficients 



3.1 The Orbit- Averaged Fokker-Planck Equation 

As described in Chapter [21 it is not advantageous to use energy E and angular momentum J 
(or its vertical compenent J^) for the integrals of the motion here as prior two-dimensional, 
spherically-symmetric studies have done [44j; any change in the gravitational potential of 
the system will cause changes in the distribution function when it is expressed as F{E, J) 
or as F{E,Jz). What is required are the radial action Ii and tangential action I2 = J, 
which are adiabatic invariants as the potential evolves - and so also is /(/i,l2) given the 
assumption of weak encounters. As given by Binney and Tremaine 02], orbit-averaging 
is then a simple averaging of quantities over the 2tt change in (an) angle variable. It is 
customary to average over the radial angle. The orbit-averaged Fokker-Planck equation 
then takes a particularly simple form: 



^/(Ii,t)--^[/C,] + ^^^^ 



/(^-^) = + o7J71Jr[/A,] (3.1) 



where f{Ii, t) = f(Ii,l2,t) is the distribution function and Cj, Dij are the drift and diffusion 
coefficients. Summation over repeated indices i and j is implied. 

Physically, the Fokker-Planck equation can be understood as a collisionless Boltzmann 
equation df/dt = with a collisional term addec|l| to account for particles (here, stars) 
scattering in and out of a given volume of phase space. When the collisional term is 
expanded in a Taylor series and truncated after the second order term, the remaining two 
terms are the drift and diffusion coefficients; as stated by Binney &: Tremaine [42j, they 
describe the expected rate of change in Ij or /j/j for a given test star with actions 
When the coefficients satisfy the relations 

Dij = Dji (3.2) 



^ "Collisional" here refers to the relatively distant gravitational encounters that dominate the overall 
dynamical scattering, and is not to be confused with the close-range "collisional mergers" of stars which are 
treated in 3331 
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and 

C, = \^D., (3.3) 
then the Fokker-Planck equation can be recast in a flux-conservative form 

(3.4) 

It win be seen below that this flux-conservative form has a distinct advantage when it comes 
to finite differencing. 

The Fokker-Planck approximation remains valid as long as number N = J dlidl2f ^ 1, 
and requires that the time step At satisfy tdyn ^ At <C ir- We have a dynamical time (as 
defined by Tremaine and Weinberg [30j) of tdyn = \/Stt/IQGp < 10"^ yr, and a relaxation 
time of tr > 10^ yr or more, and so in practice finding an appropriate size of timestep At 
is not difficult H 

Generalizing to the case of multiple distribution functions fq and adding ad hoc terms Lq, Gq 
to account for losses and gains due to stellar mergers, as well as analogous parameters Bq 
and Rq for stellar births and deaths ("remnants") due to stellar evolution, the most-general 
Fokker-Planck equation in this study is 



9 „ _ 1 9 _ d 



Qim^^) = 2dri^^J qTU\ -Lq + Gq-Bq + Rq (3.5) 



in which subscript q refers to a given subpopulation of stars labeled "g", usually distin- 
guished by being of mass ruq. 

The remainder of this chapter develops the physics of the various terms on the right hand 
side of (13. 5p . First the general perturbation potential is converted into action space for 
use in the diffusion coefficients Dij. The specific case of field-star perturbers is developed, 
which is then used to derive expressions for Dij. The other relevent perturbation potential, 
that of a stellar bar, is then constructed and its diffusion coefficients also developed. How 
to finite-difference the Fokker-Planck equation using these diffusion coefficients follows, and 
expressions for the stellar merger loss and gain terms Lq and Gq are then derived. The 
chapter closes with a discussion of the possible effects of binary heating on the model 
and the reasons for not including it in the calculation; despite its omission from the Fokker- 
Planck model, the amount of possible binary heating that could have occured is still tracked 
througout the simulations in order to confirm that the reasons for omitting it remain valid. 



3.2 The Perturbing Potential 
3.2.1 General form of the expansion 

The drift and diffusion coefficients describe the rate of transfer of action (momentum) to a 
star due to interactions with other individual stars or with a bulk stellar bar perturbation. 

^For a given test star, the relaxation time is given by tr — maxi ( -^j^ ) . In Chapter [5] tiie median 

relaxation time [7] of trq = af/ [471^/3/20^ mq\n{OAN)] is used instead, in order to indicate the relaxation 
time of stars of a particular stellar mass in position space, instead of in action-space as considered here. 
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To do so, first it is necessary to express the perturbing potential in action space. This has 
been done by Tremaine and Weinberg [30]; the following summarizes their results. Starting 
in traditional spherical coordinates when expanded in spherical harmonics is 



oo I 



^* = E E '^Ur)YU9,0)e^"^^^-''''\ (3.6) 



1=0 m=~l 



(Label "*" can refer to any perturbing mass, i.e., either to a stellar bar as a whole or to a 
single field star.) Expressed in action space, p.6p becomes 



oo oo 



^*=Y1 E ^fcnm(/l,/2)e*(*^"'l+""'2+™"'^-^**) (37) 



m=0 fc,n=— oo 

where 



°° / 2 \ 



/=0 

and in which /3 is the perturber's orbital inclination (thus /3 = for a bar in the rotational 
plane of the system). Quantity Vinmif^) contains the effect of rotating the frame of refer- 
ence so that it is aligned with the orbital plane, and Wkinmi^i, ^2) is the strength of the 
interaction of the (A;, n, m) resonance for multipole expansion term I, as will be seen belo'vJl. 
In particular, Wkinm has the form 

1 r 

Wkinmih,h) = - dwi cos[kwi - n{ip - W2)]^im{r) (3.9) 
71" Jo 

in which one may note that angles wi and — W2) are functions of r as given by (j2.6p 
and ()2.7p . By inverting (j2.6p one obtains r{wi) and by extension {ip — 1^2) {wi) as required 
for use in (|3.9p . Various symmetries of spherical harmonics and in the expansion are of 
relevance: 

• <J>* is real, and all perturbers studied have real for all {I, m), so ^i-m = {—)"^^im] 

• the system is symmetric about the z = plane, so (/ + m) must be even; 

• VinmW) = unless {I + n) is even; 

• the above two items result in Vinm being real for all (Inm), and so ^knm is real for all 
values of {knm); 

• Vi^n-m = {-)"'Vinm and W^kl-n-m = {-)"'WMnm, and SO ^k-n-m = *fcnm- 



Note that the choice of radial basis functions ^imif) is not unique [23]; the current method 
is most convenient for dealing with inclined orbits while still allowing development of the 
diffusion coefficients in terms of orbital resonances, below. 

■^For completeness' sake, we note here that Vl„m(/3) = (— )'™~"-'''^r;„m(/3)yin(-|, 0), where r;„m is the 

Slater rotation matrix = (^^gg^gl^ tan2*+™-"(|) cos^'Cf ) f3^. The sum is 

taken over all t such that the factorials' arguments are all positive semidefinite, and similarly rinm = if 
any of (l ± n) or {I ± m) is negative. 
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3.2.2 Effect of Orbital Inclination 



Many of the physical quantities that determine the diffusion coefficients depend either on 
orbital inclination /3 and/or the polar angle of the curent position 9. For the former depen- 
dence we invoke the orbital inclination distribution h{f3) developed in ^2.4.1[ Quantities 
such as Vinm{fi) are then averaged over (0,/3max) with weighting function notationally 
this will be denoted by angle brackets with subscript /i(/3). Thus for any quantity y 

^ r /3max 

{y)KP)^^ d/3hi/3)yi/3). (3.10) 

Pmax JO 

Normally /3max = tt, which accounts for both prograde and retrograde orbits. As the (3 
dependence of VinmW) is strictly trigonometric, if the form of h{/3) allows it and if no part 
of as contained in Wkinm depends on /3, the averaging of ()3.8p can be done analytically; 
else it must be done numerically. This will be of help later, when weak approximations will 
be all that are needed to satisfy the requirements for doing analytical averaging. 

A further example is given by the square of (j3.8p : 




(3.11) 



The above ^^kmn is what will be used in the diffusion coefficients below. To calculate the 
average, the squared sum is expanded and evaluated term by term; this is more efficient 
than averaging over the entire squared sum at once. For convenience, the triple (knm) will 
be written as (^i^2^3) in most sums. 

The direct calculation of ()3.1ip takes an unacceptably long time, largely due to the com- 
plexity of the calculation of ^im{r) as will be seen below; it is sped up by evaluating the 
integrands at preset values of 9 and wi , and interpolating from those values when perform- 
ing the integrations. Similarly, solving ()2.6I) for r{wi) and then integrating (j2.7p for W2 — tp 
each time they are needed in ()3.9p would be computationally prohibitive, and so they are 
also calculated on a grid and interpolated for. 

The ability to transform away the 9 dependence of Yijn{9,{)) in going from p.6p to ()3.7p 
notwithstanding, other quantities that depend not on (3 but instead on 9 directly require a 
slightly more subtle treatment. The approach here is to average any such quantity over its 
allowed range of^— /3<0<^ before doing the averaging over (3 itself. Allowing there to 
also be a purely /^-dependent factor y: 

{y; x)g = f^^^ dp h(J3)y{f3) - / ' # x(0(V', /3)). (3.12) 

Pmax Jo ^ Jo 

The inner averaging is actually performed over ip instead of 9 because, averaged over all 
stellar orbits of a given {I, J), dijj/dt is a constant whereas d9/dt is not. (See (12.71) and 
note that the phases of wi and W2 are independent; Figure [27T] gives a visual depiction.) To 
convert between 9 and ^p, (j2.8p is used. 
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3.3 Field Star Perturbers 



For an individual star of mass (not to be confused with spherical harmonic index m) at 
position (r*, 0*, cp^^ = 17=,,^ + (po) in the field, the potential at {r,9,4>) is 

$.(r;r,) = (3.13) 
|r — r*| 

oo ; . I 



1=0 ^> y ^ ) rn=-l 
oo I , I 

^ r'+' (21 + 1) 



m,=—l 

with the traditional r< and r> respectively being the smaller and larger of r and r* . Initial 
phase 0o can be absorbed into the initial phase of canonical angle variable w^, and so has 
already been dealt with in the previous section. Comparison with (j3.6p thus yields 

$zm(r;r,) = -Gm,^——-{-y-YUe,,d). (3.14) 



r 



> 



(2/ + 1 



The field star distribution is defined as the overall cluster distribution, excluding those stars 
specifically assigned to a stellar bar (which is labeled "i?" to denote that it is a bulk per- 
turber, different from individual stars) and so is expressed in terms of the canonical actions: 
/*(/(, I2) = Ylq fqi^ii^i) ~ /b(-^I) -^2)- (Here, primes indicate actions of the perturbing field 
star, and not of the stellar orbit being perturbed.) 

To deal with the undesired perturber radial coordinate r^,, we perform an additional orbit- 
averaging of the field star's contribution to the overall perturbing potential as given by 
p.l4p . this time averaging over the perturbing field star's orbital range. It is more straight- 
forward to directly average over weighted by l/|fr*| than to convert to the corresponding 
canonical angle w'l first, which would accomplish the identical task but require an additional 
step. Both Yim{Oi,.,Q) and Vr* depend on 0^, (the latter via <I>(r*)) which can at the same 
time be averaged over using the technique of ()3.12p : 

l;(ci>,„(r;r.)) = /l; ^ dr.^^ 11 ) (3.15) 




(2^ + 1) \ Jr.. [2(i5;,-c^(r,))-/f/r2] 




Of course for orbits that are very nearly circular (rp^ ~ ?'a*); ^* is well-defined for any given 
value of I[ (noting ^ 0) and so the above averaging over r^, is unnecessary; one can then 
simply take ^im{^) = ^imi''"', '^*)- Otherwise the full form of (|3.15p is required, and is what 
is used for the ^imii") factor in (j3.9p . 
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3.4 The Diffusion Coefficients 



With knowledge of the perturbation coefficients '^^^i^i^ we can now proceed to compute 
the drift and diffusion coefficients Cj and Dij. The former were derived by Van Vleck [66] 
in a quantum mechanical context, although a more related calculation is that of Tremaine 
and Weinberg [30] who derived the dynamical friction exerted on an external satellite by 
a stellar system. Our situation is analogous, but requires the generalization of considering 
both components of the action instead of only 12- Thus the derivation of Dij here follows 
a similar form. 



Denoting first-order approximations by Ai, Hamilton's equations for the system are 



(3.16) 



with generating function 



Xi 



-Re 



h=0 ei,i2=-oo 



(3.17) 



in which J7p = and uj = IsQ^+ir] is the frequency of the perturbation potential including 
a slow "turning on" of the perturber in the distant past effected by r/ > 0. Throughout 
this section repeated index p will indicate summation over all possible values 1,2,3 of the 
index, and the same character p when used in different factors indicates distinct sums, e.g., 
{epnp){£'pnp) = (£iJ7i +£2f^2 +4S^3)(A^^i +^^^2 +^35^3)- So, to second order 



Alii All,- = 



dXi dxi 



dwi dwj 



(3.18) 



Re 



h=0 li,l2=-oo 



J{ipWp~U}t) 



i{£pQp — uj) 



Re 



1^ 1^ '^j'^hMsJTJi 



Kt'pi^p - uj' 



f] sm[ipWp 



£'3=0 eJii'2=-oo 
{tpVLp — ^3$!*) cos{£pWp — £30**) 



{£pnp - £^n,Y + r/2 
7]sin{£'pWp - f^n^t) - {£'pnp - £'^n^) cos{£'pWp - f^Qj) 



where for notational convenience, a prime on a quantity indicates that any indices it takes 



32 



are primed, i.e., oj' = {i'^Q.^, + ir]) etc.. To second order, the rate of change of (A/j A/j) is 



-[Ai/i Ai/,-] = 27] [Ai/i Ai/,] + e^"* 



(3.19) 



,r/sin(^>p - fgf^.t) - (^;i7p - i'^n,)cosii'pWp - i'^nj) 



+ e 



2rit 



(i'^fip - e^n,)^ + 7]^ 

rism{£pWp — Is^l^^t) — {IpQp — h^*) cos{£pWp — l^^l^^t) 



V Sr'\-'^'ir,'(o'n ./ o N ^ cos(^;tz;p - ^l^.t) + {1'^% - ^0,) sin(^;u;p - ^O.t) 
in which for brevity ^' = ^^i^j^a ^'^'^ ^' = ^'l\l2h' 

The diffusion coefficient for an individual perturber is ^(Ai/j Ai/j), averaged over all initial 
phases of the Wp] as these are assumed to be distributed randomly, most of the terms in 
(I3.19P vanish or cancel out. Considering each of the three "grand sums" [i.e., each quantity 
contained in outer sets of brackets in l3.19p separately, this happens by the following reasons: 



1. all terms without either Ip = l'^ for all p, or else ip 
averaging; 



t' for all p, vanish upon the 



2. all cross terms of the form [±s\ii{lpWp — ^3VL^,t) cos{£pWp — £sQ,^t)] also vanish upon 
averaging over (any) initial phase angle; 

3. ip = —i'p requires that ^3 = 0; thus in each grand sum, for each term with a given 
ii = —£[ and £2 = —£'2 there is an equal-magntitude, opposite signed term with 
£1 = £'i and £2 = £'2, and so all terms with £3 = cancel out; 

4. in the second and third grand sums, the remaining terms can all be grouped into 
compound terms having angle dependence of the form [cos'^{£pWp—£3Q^,t) — sm'^{£pWp — 
£3^*4)], which also vanishes upon averaging over initial phase. 

With the above, the only remaining terms are those from the first grand sum having £p = £'p 
and £3 > 0, and excluding cross-terms. Using angle brackets to denote averaging over initial 
phases, the reduced form of the averaged perturbation is thus 



d 



(3.20) 



00 00 



E E ^^^j'^^^i.is 

^3 = l£l,^2 = -0O 



ris'm{£pWp — £s^l^t) — {£pVLp — £^0.^) cos{£pWp — £3^^,t) 



-i2\ 



.pQp - ^3J]*)2 + r/2 



1 



T] e 



2r)t 



EE 



i£pnp - £3^,)^ + 7]^ 
■2 
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in which {sm{£x + b) cos{£x + b))^ = and (^sin'^{£x + 6)) = (^cos^{£x + &)) = ^ have been 
used to ehminate the trigonomic factors. 



The analogous derivation of the drift coefficients performed by Tremaine and Weinberg [30 
found that the effect of the first-order perturbation vanishes, and that to lowest order 



-^2 



From (|3.20p and (|3.2ip it can be seen that (|3.2p and (j3.3p do indeed hold; from here on the 
flux-conservative form of the Fokker-Planck equation (|3.4p will be assumed and no further 
explicit consideration of the drift coefficients is required. 

In order to be used in the Fokker-Planck equation, the limit ?? — )• must be taken, which 
allows us to use the relation 6{x) = ^ lim^_!.o r]\x — iry|~^. Thus, one obtains 

, \ oo oo 

-[Ai/,Ai/,]\ =7r^ Yl (dj^i^i2iA^P^p-W- (3.22) 

To obtain the full diffusion coefficient one must integrate over the full phase space of the 
perturber. With /*(/(, denoting the distribution function of the perturbers under con- 
sideration and invoking the /^-averaged expansion coefficients of the perturbing potential 
from (13.1ip one obtains 



1 d2/'/^(/;,/^)^|[Ai/,;Ai/,]^ (3.23) 

„ oo oo 



^3=l^l,^2=-oo 



Note that in this formalism the effect of the dynamical friction which underlies the interac- 
tions is expressed in terms of resonances between the orbital frequencies of any star and the 
secular frequency of the perturbing potential (which is either that of the bar, or the sum 
of all potentials of field stars of a given orbital frequency ri*). In practice, for a field star 
perturber, is evaluated by using the delta- function that describes the resonance to elim- 
inate one of the action integrals and then integrating over the other. For the perturbation 
due to a stellar bar (for which fi^, is taken to be a common value for all stars comprising 
the bar), the resonance delta- function does not integrate out in (j3.23p and will have to be 
handled separately as will be seen below. 



3.4.1 Mass Segregation 

Because '^i^n^g,^ scales linearly with m^,, and is the number- weighted distribution function, 
one can see that Dij also scales with m*; this results in the same mass segregation effect as 
described by Quinlan and Shapiro [7J in which the tendency towards equipartion of energy 

^The expression here differs from the one derived by Tremaine and Weinberg [5D] who additionally had 
an integral over dj^, but in our case that's already been taken care of through the averaging over /3. Also, 
/, here is already "orbit averaged" and no summing over the ranges of w\ and is needed. 



34 



causes lighter stars to receive kinetic energy from heavier stars, which then preferentially 
sink towards the center. 

Mass segregation is observed in globular clusters, and can be qualitatively described using 
energy equipartition. The segregation can be unstable in that a population of more-massive 
stars decouples dynamically from the lower-mass population if no equipartition is possible; 
whether this happens depends on m'^/m^, and the ratio of total masses of the two populations 
|67j . Mass segregation is one example of the complex dynamics that can arise when there 
is a spectrum of masses in the system. 

As a more specific example, McMillan et al. [68] calculate that if < 2OM0, then a 
subpopulation forms in a globular cluster core in ~ 0.2 of the half-mass relaxation time frh; 
and the central density of the larger-mass stars increases by 2-3 orders of magnitude for 
"typical IMFs"; in denser clusters it is possible for this to occur before many stars evolve 
to supernovae. Even a modest range of masses can produce enough mass segregation to 
greatly increases the rate of core collapse |54j ) . The minimum central density determined 
by McMillan et al. as leading to mass-segregation-driven core collapse in globular clusters 
is similar to that used in some of the cases studied here. 

Is Mass Segregation Primordial, or Dynamic? 

Whether mass segregation in stellar clusters is primordial or is the result of dynamical 
evolution is a matter of some debate: N-body simulations of globular cluster dissolution 
performed by Baumgarft et al. [69] on the specialized GRAPE6 computing engine matched 
observed properties of globular clusters when an initial mass segregation was assumed and 
did not not match when no initial segretation was included; they also deduced a "near- 
universal" mass function for low-metalicity star formation environments. However, using 
the same "NB0DY4" code McMillan and Portegies Zwart [55] found no firm evidence for 
a priori mass segregation in young dense clusters. In Fokker-Planck (nonrotating) models 
of individual clusters, Amaro-Seoane [4_lJ found no evidence of initial mass segregation 
while acknowledging that other models for young clusters had predicted it. While none of 
the above studies looked at exact analogues of the early-universe stellar clusters studied 
here, the conservative approach is still not to assume any initial mass segregation, so that 
any mass segregation observed in the simulations is the result of the dynamical friction. 
Collisional mergers of large stars can also be a source of effective mass segregation, in that 
mergers should preferentially occur in the more-dense central region of the system. 

3.5 The Stellar Bar: General Considerations 

As stated in the Chapter [H although some angular momentum is expected to be transported 
outwards via the effects of shear between the higher- inner regions and the lower-Q2 outer 
regions, a bulk perturbation in the potential is predicted to be much more effective than 
individual stars at transporting angular momentum, as indicated by the appearance of 
the square of the perturbation coefficient ^'^^^^2^3 ™ the overall Fokker-Planck diffusion 
coefficient ()3.23p . 

Most work on stellar bars, both observed and modeled, has been on the scale of entire 
galaxies. But "inner bars" are common in observed samples of galaxies even though the 
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formation and dynamics of nested bars is still poorly understood [28]. Even for galactic- 
scale bars few studies address the issue of "strength" of the bar, in terms of the overall 
mass associated with the bar perturbation. Typical values for the strength are in the range 
of 5-15% of the overall mass available in the relevant volume: for their N-body simulations 
Weinberg and Katz [33j assumed the bar comprised 15% of the disk and halo mass within 
the bar's radius (which they took to be half the corotation radius, although their results 
were not sensitive to this choice); and Athanassoula [32] found in N-body simulations that a 
stellar bar quickly grew to ~ 5% of the mass within the radius of their stellar disk, regardless 
of whether the potential within that radius was dominated by the disk or by the halo. 

In terms of orbital composition, bars are elongated along constituent orbits' long axes |70j . 
and chaotic orbits contribute to the bar's overall makeup [32]. Chaotic orbits are not 
considered here, but this is in line with stellar dynamics as practiced since the 1940s in 
which the chaotic effects of N-body interactions are ignored in general |71jFI Preliminary 
3-dimensional modeling of stellar bars by Athanassoula |32j shows orbits similar to, but 
more general than, those of two-dimensional models. He finds that a halo component in 
the galaxy's overall gravitational potential does not stabilize the system against forming 
a stellar bar (c/. also and infers that (as yet undetermined) bar models with non- 

isotropic distribution functions would be expected to lead to even stronger bar growth. 

Dynamically, a "bar perturbation" has typically been taken as being synonymous with the 
presence of an additional quadrupole (m = 2) component in the potential {e.g., Binney 
and Tremaine [12]), and it will be taken as so here. It is worth noting that Athanassoula 
|32j has recently found some bars to have relatively strong m = 4, and in some cases even 
m = 6 or m = 8, components. These modes represent a possible avenue for further study 
but to lowest order would be expected to simply enhance the rate of transport of angular 
momentum within the cluster. 

3.6 The Stellar Bar: Implementation 

Early in the study of bar dynamics, a high central density p(0) was argued to preclude 
bar formation due to the disruptive presence of an inner Linblad resonance (ILR), and the 
"Toomre" Trot/lW^I ^0.14 instability criterion [73] was formulated as a test for whether 
or not such a non-axisymmetric perturbation could form in a system. However, high-/3(0) 
barred galaxies with ILRs are observed, and in general bars are observed to be much more 
ubiquitous than the Toomre criterion allows. Thus the true trigger of bar formation remains 
unknown [71] . Of particular interest is the finding of Gadotti and de Souza [75] that bars 
can form in Plummer sphere potentials which are embedded in nonspherical halos but which 
possess no disk component. Many studies {e.g., Athanassoula [32]) note that bars are very 
common features. 

That being said, all simulations performed here possess Trot/I^l ^ 0.17 and thus satisfy not 
only the Toomre criterion but also the Tfot/lW^I value of 0.171 calculated by Christodoulou 
et al. |76j (which was for stellar Maclaurin spheroids, but was the strictest bar-forming 
criterion found in the literature). This is not by design, but is a function of using an 

^Gurzadyan [71] goes on to state that it would be expected from ergodic theory that N-body interactions 
and the ensuing chaotic effects would dominate over 2-body interactions - however, it is emphasized that 
this applies for the case containing a central massive object, unlike the situation here. 
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observationally-motivated value for the overall rotation parameter of A ~ 0.05 as discussed 
in §2.6.21 Thus the bar is simply assumed to exist. 

N-body studies show that bars form over 1-2 rotational periods but last for many more 
10) [77], and that the stars trapped in the bar do not undergo much individual dynam- 
ical evolution [78]. Some N-body simulations show the bar losing strength over time, but 
this effect has been found not to be due to angular-momentum transfer weakening the bar 
directly but is instead due to a vertical-buckling instability of the bar in the presence of 
background gas [79], which does not pertain here. This has led us to build the nonaxisym- 
metric potential from a fraction of the lowest-mass component of the stellar population, and 
to force those stars to orbit in "lock step", sharing a common orbital frequency and phase 
(i.e., it is assumed that they all orbit at a frequency whenever the diffusion coefficient 
calculation requires knowledge of the orbital frequency). We assume there is no transfer 
of stars from bar to field or vice versa, consistent with studies showing the same stars re- 
main trapped within the bar (|80]. [78|). N-body studies show that the bar frequency is a 
compromise of the orbital frequencies of its component stars [70] , so our fie is set by either 
conserving the total angular momentum of those stars, or by simply averaging over their 
individual orbital frequencies; each possibility is dicussed below. 

Informed by the considerations of §3.5^ the most straightforward way to incorporate a stellar 
bar into our model is to assume that a fraction of the stars are trapped in the bar from 
the outset, that the resulting nonaxial perturbation in the potential is dominated by the 
quadrupole term, and that the structure rotates in bulk at a common angular speed. With 
these assumptions the stars comprising the bar can be considered as simply a subset of the 
overall stellar population for the purposes of calculating updated values of the distribution 
function / or stellar density />, while for determining the bar's contribution to the diffusion 
coefficients Dij the bar is treated as a semi-solid object with a single overall rotation speed 
JIb- Allowing the bar distribution to evolve like the field stars avoids any problem of 
inserting the bar binding energy by hand. 

The mass fraction assigned to the bar here is typically 1% of the total system mass, for 
simplicity taken from the lowest stellar mass in the IMF being used (although the choice 
of stellar mass in the bar was found to have little or no effect on the simulation's results). 
Using the definition of bar strength suggested by Weinberg and Katz [33j as the ratio of the 
bar's mass to the total mass within half of the corotation radius, this gives a bar strength 
of 5-10% depending on the initial cluster model being considered, well within the 5-15% 
range of previous studies. 

One possible concern is that the general formalism of resonant interactions represented by 
the diffusion coefficients ()3.23p is only valid when the bar's rotation speed Qb does not 
change too slowly [30j; presumably this is because if dO^^/dt is too small, the back-reactions 
of the resonant interaction on the bar itself will not be smoothly distributed over frequency 
space and may disrupt the stellar orbits comprising the bar. In this study a typical value 
for the timescale of the rate of change of the bar frequency / {dQ^ / dt) is > 10 — 200 bar 
rotation times; this timescale is consistent with the range of Weinberg's [34j calculation for 
a general bar interacting with a halo and with Athanassoula's [72] N-body bar simulations. 

The back-reaction issue may explain why, when a bar mass-fraction of > 2% (correspond- 
ing to a bar strength much greater than 10%) was tested, the simulations were unstable 
to numerical divergences in which the results were not consistent for different choices of 
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numerical stepsize At. Backing off to a bar fraction of 1% of the total system mass allowed 
for numerical stability, at the cost of being a slightly conservative choice for the strength of 
the bar perturbation. 



3.6.1 Bar Speed Determined by Angular Momentum Conservation 

Being made up of point masses, the bar's effective moment of inertia I^fi about the polar 
axis can be expressed as simply 

I^g = J dV sin^ 6 (3.24) 

in which p^ represents the smoothed mass distribution of stars in the bar. (Throughout 
this study, all quantities p are smoothed stellar densities.) The factor p^ allows for ad hoc 
tuning of the moment of intertia, and is intended to account for the fact that the bar is 
not a solid object but is a collection of individual stars streaming through the bar's overall 
structure; Pb can even in principle be negative but in general is of order Pb ^ Pb/p [ST] . 
and so normally we will set it as an input parameter. (Note that Athanassoula [32] and 
others treated the bar as a rigid object, but they did not "build" their bar models out of 
the consituent stars of the system.) 

Quantities depending on the radial coordinate, e.g., pa and e, are stored in terms of ellip- 
soidal coordinate and so it is convenient to peform the integration over a instead of r. 
The mass element is then pBdV = 4:7rpBa'^{l — e^)"*^/^ da, and one achieves 



2ttpb j dcosO j dapBa^{l-e^ f/'^a'^(^m^e + ^^-^s\Ti^e (3.25) 
Attpb j daa'^PBil - e^ f/"^ j dcosO- ^ ^ 



27r/iB / da^ a"^ Pb 



l-e2(l-cos2 0) 

1 — / sin~^ e 



e2 



e 



The above uses standard tables and trigonometric identities to do the ang ular integratioiJ!] 
[52] . The quantity enclosed in square brackets in (j3.25p approaches a value of | in the limit 
of e — 7- as is expected when finding the moment of inertia of a spherical shell; we set it to 
be exactly | for tiny values of e (< 0.001, which gives an accuracy in Xcfj of better than one 
part in 10^) in order to avoid any numerical problems of dividing by a vanishlingly small 
denominator. 

The bar frequency is then 

nB = ^ (3.26) 

where Jb is the total angular momentum of the stars in the bar, Jb = J dl dJ fsil , J)J ■ 
Note that the integral for Jb should properly be of Jz and not of the full J. We do not 
invoke an isotropy argument to approximate Jz/J because it is more straightforward to 
simply modify Pb accordingly, increasing it by a factor of ~ 2 or possibly ~ -v/27r over the 



^Substitution z = 1 — e^(l — cos^ 9) aids in the evaluation of the term with cos^ 6 in the numerator of 
the second hue in (13.251) . 
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value li^c:^ p^l p given above. 



3.6.2 Bar Speed Determined by Angular Prenquency Conservation 

Alternatively, one may assume that the bar pattern speed is simply the average over all 
stars comprising the bar of those star's orbital frequencies, as has been found to be the case 
in N-body simulations of barred systems [70]. This case is very straightforward: 

^B = ^ J (flm^UVt2 (3.27) 

As stated above, for simplicity in this study the bar is assumed to consist only of stars 
of a single mass, i.e., = mi. The total mass in the bar is simply Mb = J d?I f^. 
This method of determining fl^ has the advantage of not requiring knowledge of the bar's 
moment of inertia and so involves no assumptions about ^b- In general we employ ()3.27p 
by default; when (j3.26p is used it will be specifically noted. 



3.6.3 Bar Perturbation 



As the physical effect of the bar is largely to add an additional quadrupole term to the 
overall potential, the bar can be modeled numerically by the following three terms in the 
expansion (j3.6p . with all other terms zero: 



• the $2±2 term that describes the quadrupole interaction; 

• the $20 term that incorporates the ellipsoidal nature of the bar [83] ; and 

• the $00 term that sets the zero point of the bar's potential. 



Only the m = terms in (|3.6p contribute to the value of $b averaged over phase angle (/>; 
because of this they are useful in determining how large the m = 2 term is, despite the 
fact that they do not directly contribute to any transfer of angular momentum themselves. 
The "bar potential" <I>b is simply the fraction of the overall potential attributable to the 
subpopulation of stars that comprise the bar (as found by using in place of p in I2.4ip . 
The requirement that the expansion sum, averaged over (p, matches the known value of the 
bar potential ^b{cl,6) implies that at a given semimajor radius value a, 



$oo(a)^oo + $2o(a)>2o(^, 0) = $3(0, < 



(3.28) 



Evaluating (j3.28p at the pole and equator 
$20 and then <l>oo gives 



and ^ = f respectively), and solving for 



and 



$2o(a) 



^00 (a) 




$B (a,0) - $B (a, 



vr 



(a,0) + ^$B (a,^ 



(3.29) 
(3.30) 
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For the strength of the quadrupole, $22 is chosen so that the more elhptical the orbits of 
the stars in the bar are, the more the $22 term contributes. By the symmetry of the bar, 
the m = ±2 terms in ()3.6p are equal, leading to 

<I>2±2(a)y2±2(|,0)=i(l-|^)c|>e(a,|) (3.31) 

i.e., circular orbits contribute nothing to the bar quadrupole, and purely radial orbits make 
for a maximally strong bar. 122 is at its largest magnitude at 6 = f , so this guarantees that 
the bar potential as expressed by this expansion remains negative at all </> while its average 
over (f) vanishes. The averages in ()3.3ip are over all bar stars whose orbits cross radial 
coordinate a. When a is greater than the outmost point flmax contained in the numeric 
radial grid, each is assumed to drop off as 1/a' (or 1/r' when the underlying potential 
is spherical symmetric) from the value at Omax: ^im{0' > Omax) = (^^^^) ^«m(a)- 

Now that the strength of the bar perturbation has been characterized, the remaining diffi- 
culty is in evaluting the integrals of ()3.23p . Note that, by incorporating <1>b as defined for the 
entire bar as a whole in the calculation of the bar's (and hence ^knm) coefficients, the 
integration over the bar's distribution function /b has effectively already been performed; 
this contrasts with the field-star case in which the (and ^knm) are defined for individual 
stars and only afterwards is the distribution function integrated over. Thus (I3.22p . with its 
"bare" (5-function, is actually the proper expression to use for the bar's contribution to the 
diffusion coefficients. 

The problem then is that, for the bar, the 5-function is then no longer a function of an 
integration variable. This can be thought of as being due to that, by construction, the 
stars comprising the bar's gravitational potential orbit in "lock step", sharing a common 
bulk orbital frequency, and so the (5-function no longer has a distribution of frequencies on 
which to act. Fortunately, the fact the calculation is being done on a numeric grid imposes 
a finite scale over which any bulk resonance effect must be "smeared" anyway - namely 
the difference Af2(/i,l2) between the nearest points on the action-space grid - and so the 
resonance effect must be widened from (5-function-width to at least this size. Thus to take 
the place of the 6 function in the bar's diffusion coefficient, we insert a Gauusian window 
function J^b : 

J-b(O) = ^e-(f^-^3nB)V2(e^3nB)^ (3.32) 

so that for the (^1,^21^3) resonance, the closer the frequency Q = £iQi{Ii, I2) + £2^2(11, h) 
for any given gridpoint (/i,/2) is to ^s^^b, the more effect of the resonance it experiences. 
For the width of the resonance we compare the bar pattern speed to the timescale of the 
bar's slowing down due to dynamical friction: 

e = nB/nl. (3.33) 

According to Tremaine and Weinberg [30 1, e as defined above is of order the fractional mass 
density of the bar and so it is set accordingly as an input parameter. The normalization 
factor A = l/(-v/27re£3r2B). With this window function we can define an effective bar 
diffusion coefficient 

00 00 

^y=^E E ^i^i^W3-^B(J^). (3.34) 

i3=ih/2=-oo 
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The above Pjj is used in place of Dij in (|3.22p : this scheme allows the bar to act as a single, 
bulk perturber of frequency f^a whose effects are felt only on the "cells" in the action space 
grid over which the line of resonances 5{iii}i + £2^2 — -^s^b) falls. The model's computer 
code issues a warning if the smearing width si^i^B is smaller than the inter-gridpoint spacing 
of Ail(/i,/2) values for any action-space grid cell in which the resonance does fall. 



3.7 Finite Differencing Scheme 
3.7.1 Numerical Stability 

Any finite differencing must be numerically stable to be useful, and so here we employ the 
Crank-Nicholson method, which averages implicit and explicit schemes and is guaranteed to 
be numerically stable for any stepsize - although the diffusion coefficients, which are very 
weakly nonlinear and are very expensive to calculate, are treated explicitly. The resulting 
sparse matrix is solveable using standard methods as shown below. 

The Fokker-Planck equation (j3.4p expressed in Crank-Nicholson form is 

^ f f{T+l) f{T)\ _ 1 ^ f n fx+ly — fxy fxy — fx-ly \ ,„ „j.n 

I ^ / n ~ fxy p. fxy — fxy-1 \ 
+ 7W I -'^22xy+ —T- IJ22xy^ —7- 

'-^'^y+y- \ '-^•^y+^y ^"^j/jz-i / 

I ^ Id fx+y+ ~ fx+y- ^ fx-y+ — fx^y. \ 
+ A7 I ^12x+y Di2x.y I 

1 / 7-1 fx+y+ ~ fx-y+ „ fx+y- — fx 

" Aj — I — a7 — a7 — 

'-^'Jyj^y— \ '-^^x+x— '-^■'-x+x— 



+ 

The notation in p.35p requires explanation: superscripts denote the current already-solved- 
for timestep (r) and the to-be-solved-for timestep {t+1); they are implicit for all / values 
on the right hand side. The averaging over (r) and (t+I), i.e., over the explicit and implicit 
differencings of (j3.4p . constitutes the Crank-Nicholson scheme. Subscripts xy denote points 
on the two-dimensional grid in action space, along the radial (first subscript) and tangential 
(second) directions. A trailing "it" indicates the position partway between the labeled 
gridpoint and the one ±1 gridpoints away {e.g., D22xy- = ^22x(j/-i)+ if "partway" is taken 
to be "halfway", as will be discussed in §3.7.3p . Thus (j3.35p is an example of differencing 
the flux-conservative diffusion equation "as it stands" [84J. Finally, AI^z' = {Iz — Iz') and 
similarly for ^Jzz'- 

Note that the diffusion coefficients are treated purely explicitly; as are p etc., they are 
only dependent on the overall ensemble / and not upon any particular f^y, and so are 
similarly considered to be part of the "snapshot" of the system at timestep r over which 
the evolution of / is calculated. 
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3.7.2 Time Splitting 

Operator splitting ^84j (also called "time splitting" when referring to a time-evolution op- 
erator such as used here) is used to evaluate each of the lines in (j3.35p individually. Using 
fractional increases in r to conceptually denote this splitting, this can be written as 

A^+^) _ At) . , / f _ f f - f 

Jxy J xy -L \ 1 j ^ J x+ly J xy J xy Jx 

'At ~ 4 

(r) 

Ar+I) _ (r+i) (r+l) f _ f f _ f 

Jxy Jxy -L ^ I J xy+1 Jxy ^ Jxy J: 



At 



— = 4 E V '^''^^^ (^12x+j/ + -D21xy+ ) - fx+y. {Dl2x+y + D21xy_ ) (3.36c) 

— fx-y+{Dl2x-y + D2lxy+) + fx-y-{Dl2x-y + -C21xj/_ ) 

in which for convenience = (A/^^a._ AJ^^^ ) has been used and terms rearranged in 
(|3.36cp . Collecting terms for (r + i) on the left and for (r) on the right, (I3.36ap becomes 

_ 1 Dllx^y ^(r+|) / AJa;^a;„ 1 ^llx-+y 1 ^ Jr+l) _ 1 Dux+y_ Jr+l), , 

4A4._/--i^ +V At +4A4+i. + 4A4._J^^^ 4 A/,.+i/-'+i^^ ^'-''^ 

^ 1 Dllx^y At) f ^Ix+x^ 1 Dux+y 1 ^llx-y \ ^(^) 1 -Pllx+y ^(r) 

4 A4 x-1 ^"-'^ ^ V At 4 A4+1 X ^AIxx-iJ -^"^ 4 A/,+i , 
which can be seen to form a tridiagonal linear set of equations in x that can be solved 

(t+-) 

for the unknowns fxy ^ , for a given stepsize At and (fixed) y. Routine TRIDAG from 
Press et al. p?] is used to solve (j3.37p for each y value in turn. Given knowledge of the 

fT-|- — ) (t I ~ ") 

full fxy ^ , the solution of ()3.36bp for fxy ^ is entirely analogous, with the radial and 
tangential coordinates swapping roles. 

Time splitting is not mandatory for setting up differencing, but it greatly simplifies the 
scheme both conceptually and mathematically. It will be shown in ^4.61 that time splitting 
does result in a more accurate and robust differencing. 



3.7.3 Ensuring a Positive-Definite Distribution 

In addition to numerical stability, when tracking the evolution of a distribution function of 
real objects it is greatly advantageous to apply a differencing scheme which guarantees a 
positive-definite solution for every gridpoint and at each timestep. For the one-dimensional 
Fokker-Planck equation there exists the Chang-Cooper spatial differencing scheme which 
(when combined with the Crank-Nicholson time-differencing described above) has both 
these qualities [85]. The Chang-Cooper scheme consists of a method to formulate a working 
prescription for the "half-grid" points such as fx+y'- for example, if one defines a parameter 
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using fx_^y = (1 — 6x(^y-^) fx+iy + 5x{y)fxy then Chang-Cooper provides a calculation for 
^x{y) which guarantees / will remain positive-definite. 

A one-dimensional differencing scheme does not necessarily generalize to two or more di- 
mensions, however: Chang-Cooper does not give a general form for calculating both 6x{y) 
and 5(^x)y Fortunately when the Fokker-Planck equation is specifically cast in its flux- 
conservative form ()3.4p . the Chang-Cooper method reduces trivially to the case of "centrally 
weighting" all derivatives that involve the distribution function: = \ s-nd 5[x)y = \ 

for all X and y. Thus fxj^y = \{fx+iy + fxy) etc. for the half-grid values of /; substituting 
these into (j3.36cp and recollecting for the various values of / on the gridpoints (i.e., fxy 
and its eight immediate neighbors) yields a difference equation for fiy^^^ analogous to that 

of (j3.37p for /iy^^''. The ensemble of such equations for all values of x and y forms a sparse 
matrix system which is solved using the routine LINBCG from Press et al. |84j . 



3.7.4 Numerical Boundary Conditions 



At the edges of the (I, J) grid the differencing schemes as described above require knowledge 
of fxy values beyond the grid proper. Denoting the highest values of x and y on the grid 
as X and Y, for x = 1 or x = X and for y = 1 or y = Y numeric boundary conditions 
must replace (|3.37p and its analogues for (r-|-|). Using the x = X case as an illustration, 
possibilities include: 

fiy''' = fx-^y (3-38a) 

/iV*' = f'x^r (3.38c) 

These are all straightforward extrapolations; cases (a) and (c) are taken directly from 
Strikwerda [86], while case (b) has been generalized for a non-uniform grid. 

An option which does not require any extrapolation consists of employing a one-sided (to- 
wards the interior) differencing on the boundary instead: 

^f{Ix) _ fXy - fx-ly , 

-AT- - -AI]^ 

which can then be used in place of Chang-Cooper's centered differencing to derive analogues 
of equations (|3.35p and (j3.36p for use on the boundaries. This scheme effectively causes 
A^//A/^ and A^//AJ^ to vanish on the boundaries, and so is referred to as the natural- 
spline method (although the full differencing of the Fokker-Planck equation only vanishes if 
e.g., Diix+y = Diix_y as well). It is not a numerical boundary condition per se, but rather 
a new differencing scheme for the boundary points that avoids the need for any numerical 
boundary condition. 

Experimentation with all four possible boundary condition schemes showed that only the 
natural-spline treatment of (I3.38dp accurately tracked an analytically-solvable test case, as 
shown in Figure 14.121 and described in ^4.61 
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3.8 Stellar Mergers 



3.8.1 Rates of Loss due to Mergers 

Modeling stellar mergers due to collisions requires the transfer of stars from one fq to 
another, e.g., if a star of mass mg with actions {I, J) and one of mass rrigi with actions 
(/', J') combine into a new star of mass ruq = rriQ + iJiq', then J) and fq'{I' , J') must 

be decreased, and fq increased at a dynamically appropriate value of the action (/, J)Il| 

This alteration of fq, fQ and fq' is accomplished via the loss (Lq) and gain (Gq) terms 
added to the Fokker-Planck equation, as in (|3.5p . Note that these terms are not part 
of the differencing scheme, as they are simple rates of change and so can be calculated 
directly. The terms Lq and Gq are similar to those developed by Quinlan and Shapiro [T], 
who used a one-dimensional distribution function F{E) only. We adapt them to the two- 
dimensional f{I,J), calculating the probability of collision given cross section a. (Context 
should prevent any notational confusion between the cross section and the components of 
the velocity dispersion, e.g., u^.) Expressing the cross section as a sum of powers of the 
collisional speed |Ai;| with mass-dependent coefficients, i.e., 

J^a„(g,(?')|At;r (3.39) 

a 

then the rate of collisions of a given "target" star of mass ruq and actions (/, J) at radius 
r with all stars of mass niq' is 

K'ir,I,J) = ^aM)\Av\ = E%^ l^^'Ml',J')\^-r^' (3.40) 

where primes are used to indicate quantities dependent on the actions being integrated over, 
and the only dependence upon the target star's actions (I, J) is in Av = Iv — ifl. The above 
makes use of the density (|2.2ip as an operator, but in doing so introduces an ambiguity due 
to the lack of full information about r, v and v' in the orbit-averaged distribution function 
fq'. This ambiguity will be dealt with below. 

Given values for the actions, the magnitudes of the radial ([^^(r)!, |f^(r")|) and tangential 
{\vt\ = J/f, Wj'l = J' /r) components of the stellar velocities are well-defined, but there 
is ambiguity in the latter's subcomponents due to lack of information regarding J^/J. A 
conservative approach is to assume that the relative orientation of the tangential parts of 
the velocity vectors is random, i.e., if 7 is the angle between vt and v'rp, so that 

\Av\l = {\vr\ T K\f + (WtI - Kr|cos7)2 + |wy|2sin2 7 (3.41) 

then 7 is randomly distributed. This is conservative because it assumes no "collimation" 
of Vt and ifrp due to the effect of overall cluster rotation. Taking the average over random 
angle 7 and over the possible relative signs of the radial components one can define for any 

''In practice there isn't usually a mass bin q such that m, — rriQ + rriqi exactly, and so interpolation over 
the two bins nearest (mg + m,/) is required; this is discussed in ii3.8.3l 
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given function (7(|Ai;|), 

{g{\Av\))^ ^l-J^jY^a (I A^1±) . (3.42) 

The r dependence is easily eliminated by integrating over the orbit of the target star. Thus 
the orbit-averaged rate of collisions for a single star is 

(3.43) 

where Mqq' is given by (I3.40p . in which the endpoints of the action integrals therein are 
those energetically allowed at r. In practice the integral over dr above is brought inside the 
one over dIdJ, in which case one integrates over all (/, J) and restricts the integration over 
r to the range, if any, over which the orbits overlap. A possible simplifying assumption for 
the integral in (j3.43p is due to the fact that the target star spends the bulk of its time at 
or near its orbital endpoints Vp = r(0) and = r(7r).) 

With this in hand the total loss rate due to stellar collision of all stars with given values of 
action (/, J) and mass niq is simply 

Lq{I, J) = fq{I, J) Nqq,{I, J). (3.44) 

9' 

Of course, in reality there is a continuum of masses and so the sum in (|3.44p should properly 
be an integral, but in this study stellar masses are defined on a finite grid of values and so 
the sum is over those values. 

3.8.2 Rates of Gain 

The calculation of the gain rate Gq{I, J) is conceptually similar but slightly more cumber- 
some, owing to the difficulty in determining the new velocity vector of a star which is the 
product of a merger of two smaller stars. In Gq, I and J now denote the actions of the 
star of mass ruq produced by the collision, and so we introduce the notation (X, J') for the 
target star's actions. Thus at any given r, conservation of momentum requires 

= [{mq — mqi)J + ruqij' cos 7]^ + m'^/J''^ sin^ 7 (3.45) 

and 

mqVr{r) = {rUq — mq')i9± ± niqiv'^ (3.46) 

in which 'd± has been used for the target star's radial velocity and the choice in relative sign 
accounts for the possibilities that the stars have like or opposite radial velocity directions at 
the time of collision. Note that this is the same "it" as in (|3.4ip . although of opposite sign, 
and does not need to be independently averaged over. After solving (j3.45p for J = r\^T\ 
and (|3.46p for ■d±, the target star's orbital energy is simply 8± = \{'&'± + |i?Tp) + ^{t), and 
then X± can be found directly from (12. ip . 

The analogue of (j3.43p in the merger-gain case is the rate of producing stars of mass rUq 
and actions (/, J) from those of masses tuq = rUq — niq' and niq' and with respective actions 
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(X, ^7) and (/', J'). Because we are "back-constructing" what X and J are from knowledge 
(/', J') and (/, J), the orbit-averaging is only over that portion of the orbit of a star with 
actions (X, J') that overlaps with the product star orbit as determined by its actions (/, J). 
(Consideration of the overlap of the orbit having actions (/', J') with that of (/, J) is already 
implicit in the use of Pq'{r) in (j3.40p ). Labeling the target star's radial frequency as Ow, 
its associated canonical angle variable as Wi, and its orbital endpoints TZp and TZa, 



max 



N'{I,J) = / dw,Ma,'{r{wi),I,J)) (3.47a) 



Qq' 



7 



= ( / -iYT^faq'ir,I,J)) . (3.47b) 

Note that in this case we integrate over dWi so that contributions to N'^^, (/, J) over the 
entire orbit of the merging stars are considered, whereas in the merger-loss case we orbit- 
averaged over dwi in order to eliminate the r dependence. 

To determine >Vi(r) (or TZp{I, J) and TZa{^, -J)) one must solve (|3.46p and (|3.45p - which is 
only possible within the integral over dWi (or over dr), and which in turn requires knowledge 
of Wi(r) in order to determine the integral's endpoints. To avoid the impasse, the form 

/ i-mm{ra,r'a) 0^^,dr ■ \ 

is used instead, with the additional requirement that only dynamically-allowed solutions of 
(I3.45P are allowed: i.e., for a given value of 7 we define AAgg' (r, X, J") = if J" < or is 
complex; this is equivalent to only considering orbital ranges that do indeed overlap. 

The rate of gain due to the mergers is thus 

G,(/, J) = ^^iV^^,(/, J)/s(X, J) (3.49) 

with an extra factor of ^ inserted to counteract the effect of double-countingH it is arbitrary 
which is called the "target" star, and it takes one of each to create a merger product. The 
right side of (j3.49p contains a notational sleight of hand: since (X, J7) depend on /', J', and 
7, the integrals from (13.40p and (13.42P now act as operators on /g(X, J"). For clarity, the 
full expression for the gain rate is 

G,(/,j) = x^^^y dij^^^^ ^>_^/,(x±,j)iA.i± 

(3.50) 

in which the final sum is over the two possible sign choices in ()3.4ip . and over (one or both) 
dynamically-allowed values of J found from solving (I3.45p . In practice it is required to do 



*Other studies, such as [7], have included a factor of 1/(1 + Sq^i) instead of 1/2 to prevent double- 
counting. The distinction is that here we sum over all possible coUisional pairs and attribute their products 
to the appropriate mq = tuq +m'q, as opposed to predetermining the product m value being considered and 
invoking a delta function S{mq — tuq — m'q) to restrict the mass values of the colliding stars as needed. The 
present method is more efficient computationally as all Gq for a given (7, J) are found simultaneously. 
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the orbit-averaging integral directly over r and not over Wi, as in ()3.47bp and (|3.48p . 



An artifact of this approach is that the symmetry between "target" and "object" stars is 
not obvious in (j3.50p . However as the numerical calculations for Lq and Gq are performed 
by independent sections of computer code, a match between overall merger losses and gains 
shows the validity of the method, as will be discussed in ^4.51 



3.8.3 Mass Bookkeeping 

The above derivation made use of the manifest fact that any star which is the product of a 
purely-inelastic stellar collision and merger will have a mass equal to the sum of the colliding 
stars' masses: niq = niQ + niqi. But the model only tracks stars of a finite number mass 
values, and so it's possible that ruq does not correspond to a stored value. To account for 
such intermediate mass values, merger products are interpolated across the nearest values 
below and above niq on the mass grid [87|. Referring to these values as mq_ and rug^ 
respectively, one has 

= ~ "^'^ Gq (3.51) 
mg_^ - mq_ 

and 

= Hh^JI^Gq. (3.52) 

ruq^ - ruq^ 

Note that because actions are conserved in collisions, no interpolation in I or J is required. 
Thus the merger-gain rates from both (j3.5ip and (j3.52p for a given mass value contribute 
to the full Fokker-Planck equation (13. 5p . 



One may ask why do we not merely round down the merger-product's mass to the next 
lowest mass bin value, in order to account for mass loss during a not-completely inelastic 
stellar collision which produced the merger. The main reasons are (1) that it would also 
inadvertently reduce the total mass of the cluster, which while it would be a reasonable 
effect in an ordinary globular cluster would not be expected to occur in these more massive 
systems; and (2) that there is no way to account for diffuse gas in the calculation of the 
cluster's gravitational potential, nor in associated dynamical quantities. Thus all collisional 
mass is assumed to go into the newly-merged star, a similar approach to that taken by 
Quinlan and Shapiro [7] who surveyed results of hydro dynamical simulations of stellar 
collisions - most of which find a maximum mass loss per collision of < 11% - and concluded 
that the average mass loss will be much less than that and thus also assumed complete 
coalescence into the produced larger-mass star. It is also consistent with the finding of 
Freitag et al. fT9] that assuming coalescence is "fully justified" for velocity dispersions 
Vrms ^ 300 km/s - a condition largely satified by the simulations here, the vast bulk of 
which have a calculated initial velocity dispersion of ai < 310 km/s. 



3.8.4 The Delta-function Approximation 

The full merger-gain coefficient calculation (13.50p is by far the most computationally inten- 
sive part of the model, and so calls for a simplifying approximation. The "delta-function 
approximation", in which all stars at a given of radial distance are assumed to share a 
common = Vrms{f)-, was found by Quinlan and Shapiro [Ij to provide large performance 
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gains with only a marginal loss in accuracy in most cases. The key to the approximation 
is to not replace the density pqi in (I3.40p with its distribution-function integral form (I2.2ip 
and instead to substitute \v'\ with Urms('') throughout, thus avoiding two computationally- 
expensive integrals over the action components. However the present two-dimensional study 
requires velocity components v[. = \v'\ cos9y = VrmsCosOy and v'j, = \v'\ sinOy = VrmsS'mO^, 
and so the approximation comes at the cost of an averaging over sin 9^ : 

(3.53) 

The actions at a given position r(Wi) are determined by the new conservation rules 

rriqV^ = [{iTiq — mq')\'dT\ + mqiv'j.^gSm9v cos7]^ -|- m^/U^^gSin^ 9^ sin^ 7 (3.54) 

and 

mqVr{r) = {niq — mqi)'d± lb m qjv'j.j^g COS 9 y (3.55) 

which, as before, may not have a solution for either or both {Xj^^J) or (X^^J) (in which 
case fQ{X,J) = 0). The notation v'j.^^{r) indicates the rms velocity specifically of the 
subpopulation of stars of mass g', making this an even less-severe approximation than 
in the one-dimensional case. Note that, consistent with the conservative assumption of 
a randomly- varying 7 in §3.8. H v' is also taken to be evenly distributed over 9^. The 
averaging is restricted to the range < ^t, < ^ because the "it" in ()3.55p already accounts 
for f < 6*^, < vr. 



3.8.5 The Cross Section 



Also following Quinlan and Shapiro [7], for the stellar mergers we take a hard-sphere cross 
section with gravitational focusing correction: 



a = 7r(r^ + ) 



1 + 



2G(m + m') 
(r* -I- ri)|Au| 



(3.56) 



in which m and m' are the individual stellar masses, and and are the physical stellar 
radii of main-sequence stars of those masses. 

In practice, either term in the above may dominate in a given collision, depending on the size 
of the relative speed of collision |AiT|; Ebisuzaki et aJ. [.14j has found observational evidence 
that the gravitational focusing term can be important in compact stellar clusters, and 
Freitag et ai. [19| claim that the gravitational focusing term dominates when the velocity 
dispersion ai < 300 km/s, which is the case for most but not all simulations presented here 
(the exception being the "E4A" models - which, as mentioned above, almost satisfy that 
criterion). 



3.9 Binary Mergers, Binary Heating 

It will be seen in Chapter [5] that two-body mergers, while in some cases fairly frequent, 
do not dominate the dynamics in our simulations. Given that 3-body collisions are much 
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more likely to form binaries than to result in a direct merger, that Quinlan and Shapiro [7] 
found that few hard 3-body binaries formed during most of their simulations, and that in 
the absence of a central massive object Preitag et al. [1^ found that 2-body mergers start 
before 3-body binaries can form, we ignore the effects of 3-body collisions for stellar mergers 
as well as for the effect of binary hardening on the overall energy budget of the system. In 
addition, following the lead of prior studies ([56], |22j ) the effect of primordial binaries is 
also not considered for simplicity's sake. For determining the rate of stellar mergers this is 
a conservative approach, as when there is no massive central object primordial binaries are 
likely to foster collisions (even if once a massive black hole does form the binaries then serve 
to "grind down" stars instead of growing them) p!9|. In their N-body simulations Portegies 
Zwart and McMillan [9j also found that 3-body encounters (binary-l-star) increase the rate 
of mergers, in contrast to prior assumptions. 

Binaries can also serve to heat the overall distribution of field stars as the binaries harden 
[7], with the potential to eventually halt and reverse core collapse. However, in practice for 
dense clusters this effect is not important prior to the late stages of core collapse and so 
can be ignored in earlier stages [22j. Also, stellar collisions between binaries and other stars 
lead to a significant reduction in the amount of heating produced [56j . Finally, binary -|- 
field star interactions are most likely for very small {6v ~ 5 km/s) relative velocities [88], 
whereas the models employed here feature velocity dispersions measured in the hundreds of 
km/s. 

For all the above reasons, binary heating has not been incorporated into the simulations 
performed. However, a test of the maximum possible amount of binary heating the system 
could have achieved was performed for a variety of models, the results of which are given 
in 
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Chapter 4 



Validity Tests and Model 
Parameter Choices 

4.1 Overview 

There are five major calculational components to the simulation method: 

1. updating the gravitational potential $(r) and stellar density Pq{r) given a new distri- 
bution function fq{I, J) at each timestep; 

2. calculating the resulting new values for stellar-dynamical quantities STich as the orbital 
frequencies Qj{I,J), as well as the conversion between orbital angle wi and radial 
position ri^wi) for a given actions (/, J) etc.); 

3. calculating the full set of diffusion coefficients Dij{I, J); 

4. calculating the rates of mass loss {Lq{I, J)) and gain {Gq{I, J)) for each stellar mass 
value niq within the distribution functions fq{I, J); and 

5. finite-differencing the Fokker-Planck equation to update each fq{I,J) using D^j, Lg 
and Gq. 

Also requiring consideration are choices for the purely-numerical parameters of each simu- 
lation: 

• the size and range of the grid of actions (/, J) on which many quantities are calculated; 

• the size and range of the radial-coordinate r grid for p{r), $(r) etc.; 

• the choice of timestep At; 

• the number of terms used in the spherical-harmonic expansion of and 

• how to "bin" the range of stellar masses studied into discrete values niq. 

Tests of each of these aspects of the overall calculation are presented in turn in this chapter, 
as is verification that the effect of binary heating is sufficiently small that it need not be 
included in the model. 
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Figure 4.1: Demonstration of the ability of the potential-solver to converge even when given 
an initial guess for $ and p that are extremely in error. The solid line is the analytical curve 
for $(r) for a Plummer sphere, the long dashes are the initially input values, and the short 
dashes are the iteratively solved-for potential. This convergence was obtained after only 4 
iterations. The units in this figure are arbitrary for testing purposes. 

4.2 Potential-calculating Tests 

The first part of the model iteratively calculates the new density p and gravitational poten- 
tial $ given an updated /; it is most easily tested by giving it an artificially bad (i.e., far 
from correct) initial guess for $(r) and and verifying that it does converge on the true 
values. That it does can be seen in Figures [4 . 1 1 and [ 4 . 2 i Of note is that the input functions 
were close to being "maximally" bad, i.e., any much larger discrepency from the correct 
values produced unphysical results in intermediate calculations {e.g., p < 0). Yet in these 
tests cases the solver converged on the proper solutions in approximately the same number 
of iterations as it does when used in the actual model calculation. 

4.3 Dynamical Tests 
4.3.1 Orbital Frequencies 

The orbital-dynamic quantities described in ^2.21 appear in almost every higher-level cal- 
culation in the model. In order to test the accuracy of the calculation of the dynamical 
quantities' values, two different potential-density pairs for which all quantities are also ana- 
lytically calculable were used as test cases: the two-dimensional Simple Harmonic Oscillater 
(SHO), and the isochrone potential |42j . 

Figure [431 shows a sample result of the SHO case; the upper curves are of Oi, lower are of 
fl2- The results do not depend on whether the were found by directly integrating over 
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Figure 4.2: The same test of the potential solver as shown in Figure IITTI but now plotting 
p(r), again in arbitrary units. The initial guess was chosen to be much smaller and flatter 
than the true solution. 

the orbit {i.e., using (j2.3p or (j2.2p ) or by calculating the partial derivative of the energy 
(using ()2.5I) ): this is consistent with the fact that the numerically-determined values for E 
and for the orbital endpoints matched the analytic values to within 10~^ or better. As 
can be seen, there is a slight systematic offset in the value found for l^i, but otherwise the 
match to the analytic curves is good; this is a typical result. 

When testing the calculation of dynamical quantities, there was a choice in what values 
to use for the potential $(r): either numerically-determined values {i.e., a "fully- numeric" 
test which employed the entire calculational machinery) , or the analytically-known potential 
(a "semi-analytic" test which isolated the dynamical calculations only). Figure [13] shows 
results for both, in the case of the isochrone model. Again there is a slight offset between 
the fully-analytic 0,i and the numeric curves. Of note is the close match between the semi- 
analytic and fully-numeric calculations; this is not surprising given how well the potential- 
finding algorithm described in Chapter [2] works, as seen in §4.21 

4.3.2 Orbital Angles 

As stated in §2.21 the diffusion coefficients do not depend on the canonical angle variables 
Wj, nor is knowledge of them needed for any of the other calculations in the model; this 
is largely due to the fact that the sinusoidal dependence of the multipole components of 
the potential in action space {i.e., as shown in (j3.7p ) is eliminated upon orbit-averaging0. 
However, r{wi) and [■0 — W2]{r) are explicitly required to perform the integral of ()3.9p in 
the diffusion coefficient calculation, and 'W2{r) is needed for many of the tests described in 
this chapter. 

^Although as shown e.g., in (|3.47|l . if desired it is always possible to change variables from r to wi by 
using (|2.6p . 
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Figure 4.3: Orbital frequencies (upper plot) and Q2 vs. J for a Simple Harmonic 
Oscillater potential, at a constant energy. Straight lines are the analytic values, varying 
ones are the numerical- found solutions. Numerical values of other dynamical quantities 
(energy, orbital endpoints) are accurate to within 10"^ or better of the analytic values. 
As explained in the text, numerical curves in this plot were found using the simulation 
code's calculations of the potential and density; substituting analytically-known potential 
and density values into the calculation produced identical curves. Similarly, whether the 0,j 
were found by integrating over the orbit (using (|2.2p or ()2.3p ) or by taking = ^ (j2.5p 
made no difference. 
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Figure 4.4: 0,i (upper plot) and ^2 vs. / for an isochrone potential, at a constant energy. 
Similar to Fig. 14.31 but here the numerically-found 0, curves show both the semi-analytic- 
and fully-numeric-potential cases. 
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Figure 4.5: Orbital angle wi{r) interpolated from values calculated using (j2.6p v. the 
analytically-known value of wi in an SHO potential is plotted as a dashed line. Results of 
a semi-analytic test were identical to the fully-numeric results shown here. The solid line 
shows the diagonal, for comparison. 

Figure 1431 shows the numerically-calculated wi{r) versus the analytically-known function 
wi{r) in an SHO potential. For the calculation of the curve in Fig. 14.51 first wi{r) was 
integrated using (j2.6p at a finite number of values of r and then the shown curve was inter- 
polated at other points from that grid of values. The technique of calculating a complicated 
function on a finite grid and then interpolating to determine the function's value at an 
arbitrary point was employed throughout this study in order to reduce computation time, 
and so this serves as a test of that procedure as well. 

Figure [461 is the analogous plot of r, i.e., in which wi{r) was calculated on a grid and then 
the inverse funtion r{wi) found by doing a reverse interpolation. In this case the inter- 
polation does more than merely reduce the required computation time, it is also required 
to actually invert the function. Figures 14.61 and 14.51 show that the results of the numeric 
calculation of wi and r from <I> and p match the true values of those orbital angles. 

The other interesting orbital angle term is — 4') ^ described in §2.21 Figure 14.71 plots 
{■W2 — ip) calculated from (j2.7p versus the analytically-known values in an SHO potential, 
for two different values of the actions (/, J) {i.e., across two different orbits). The orbit 
plotted with "+" in that figure is a typical result, while the one shown with "x" is one of 
the more extreme cases. It can be seen that although (^2 — V') isn't calculated quite as 
accurately as wi (perhaps because of the delicacy of the [^2 — J/r^] factor in (j2.7p ) neither 
orbit's values are far from the expected diagonal. 
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Figure 4.6: As in Fig. 14.51 but here the interpolated inverse function r{wi) is plotted v. 
the analytic values of r. 
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Figure 4.7: Angle term {w2 — ip){r) plotted as numerically-calculated vs. analytically-known 
values. Points denoted by are for an orbit with intermediate values of actions (/, J) 
for the SHO potential used, while "x" are for an orbit with {I, J) near the upper end of 
the range of actions that are still bound in the potential. This figure is done as a scatter 
plot because, unlike wi{r), {w2 — i^)ir) is not a monotonic function of r. As in Figs. 14.51 
and 14. 6[ {w2 — tp) is first integrated from (|2.7p for a finite number of r points, and then the 
values plotted here interpolated are from that grid in r. 
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4.4 Diffusion Coefficient Tests 



4.4.1 Low-level Calculations 

The field-star-orbit averaging of the potential's expansion terms as given by (j3.15p is one of 
the more complex steps required in calculating the diffusion coefficients. But for some values 
of the expansion index I, and over some ranges of r/Rc, the integral in (|3.15p is analytically 
solveable in the test potentials - specifically, for Z = or 2 in the isochrone potential, and 
for < Z < 6 in the SHO potential. In both cases a consistent good match, within 10% 
or better, was found between the numerically-integrated values and the analytic one, and 
a match within 2% or better between numeric integrations using analytic or fully-numeric 
integrands. (Here "fully- numeric" means using only code from the model itself, with the 
test potential values as inputs.) While encouraging, this test only verifies that the integrals 
in ()3.15p are being performed accurately, and not that they are correct in their form. For 
that, see below. 

4.4.2 Diffusion Coefficients: Reproducing the Potential 

As shown in ()3.23p the major part of the diffusion coefficient calculation consists of deter- 
mining the expansion coefficients ^i^i^i^ of the potential in action space, whose relation to 
the (position-space) potential of a given perturber is given by (|3.7p . Thus a fundamental 
test of the validity of the diffusion coefficients is to check whether they can be used to recon- 
struct the overall gravitational potential of a field star or of the bar, effectively reversing the 
change of coordinates from position-space to action space from which they were calculated. 
To reconstruct the potential, first we have to add in the I = and / = 1 terms to (j3.7p that 
don't contribute to the dynamical friction: 

oo oo 

^W = Z1 '^e,i2i3il^ J) cos{ipWp-ujt). (4.1) 

'3=0 li,l2 = — 00 

Not all of the above terms are needed: as stated in ^2.21 the angle contains a random 
phase (j)o and so only terms with Z3 = contribute to the sum, due to the cosine factor. 
Angle W2 includes ip, which is measured from the ascending node to the current orbital 
position and is unconstrained by where in the allowed range of orbital radius r the star 
currently is; however, this does not imply that terms with /2 7^ do not contribute: even 
though -0 and wi are uncorrelated, '0 and W2 may be and so W2 does not act as a random 
phase as W3 does. 

So upon orbit-averaging, the cosine factor restricts contributions to only those terms with 
I3 = 0. Note that the cosine factor does not appear in the actual dynamical friction 
calculation, and so non-zero I3 terms do contribute there - in fact, their presence is what 
allows for the resonant interaction between field star and perturber. All values of li and I2 
must be included here however, as wi, W2 and the ^ are all dependent upon r for a given 
value of the orbital action vector (I, J), and so the only terms that contribute to ()4.ip are 

00 

$(r) = "^hhoil, J) cos(/iu;i + I2W2 - cot) (4.2) 

ll ,l2 =—00 
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Figure 4.8: Scatter plot showing the results of using (|4.2p to reconstruct the potential 
<I>(r) = l/(|r — r*|) of a single star within a larger cluster from its diffusion coefficients 
^hho- The various points plotted represent different choices both of test-point location r 
and of object-star action values (/, J). Units are arbitrary but equivalent along both axes. 

which can then be averaged over any object-star orbits being considered. 

The simplest test case to consider is that of a single object star's potential as given by 
$(r) = l/|r — r*|. Figure shows the values of <l>(r) reconstructed using (14. 2 p from the 
diffusion coefficients of object stars with various orbital actions (/, J), calculated at several 
locations r within the overal cluster. In general the match between the reconstructed $(r) 
and the expected value l/{\r — r^\) is fairly close, with some exceptions at larger values of 
l/(|r — r^\), i.e., when the gravitational effect of the object star is greatest. This is not a 
large cause of concern however: Figure shows that the majority of the outliers in Fig. 
14.81 are due to object stars on nearly-circular (/ ~ 0) or nearly-radial ( J ~ 0) orbits; as 
such orbits are nearly unpopulated in the actual distribution functions /(/, J) used in the 
simulations they do not contribute strongly to the overall diffusion coefhcient values. For 
the non-extreme orbits that make up the bulk of /(/, J), the diffusion coefficients are seen 
to be able to reproduce the perturbing potential as expected. 

4.4.3 The Bar Perturbation 

The calculation of the bar's diffusion coefficients, as shown in ()3.34p . was constructed so 
that the bar potential should describe a quadrupole while still averaging to the underlying 
aggregate potential of its constituent stars. Figure 14.101 shows that this calculation works 
as designed: near the cluster's equatorial plane the bar potential, as reconstructed from the 
diffusion coefficients using (j4.ip . is dominated by the $22 term, but towards the pole the 
9 dependence of the spherical harmonic I22 allows ^>oo and ^>2o to contribute most to the 
overall bar potential. 
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Figure 4.10: The gravitational potential of the bar v. azimuthal angle (/>, as reproduced by 
summing (j3.7p using the diffusion coefficients of (j3.34p . This plot is for r = 1.079pc using 
the "ElB" model. The flat line is the base axisymmetric potential that the same stars would 
have if they were not assumed to be in the bar. The curved lines show the bar potential at 
different values of polar angle ^, from Q = 0.002 for the least-curved line through Q = 7r/4 
and = 7r/2 for the most-curved. 
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4.5 Merger Losses &6 Gains 



The calculations of the rates of mass loss and gain due to stellar mergers, Lg (as given by 
(j3.44p ) and Gg (by (j3.50p ). are done independently, offering the opportunity to compare 
them for consistency. Figure 14.111 shows the results for the three models primarily studied: 
the E4B model with Kroupa- or Arches-style IMF, and the E2B model with Arches-style 
IMF. The results shown are taken from the same simulations described in detail in Chapter 
[5] and so have upper mass limits of 125Mq and either 9 mass bins (for the Arches-style IMF 
given by (12.55P ) or 10 bins (for the Kroupa IMF of (I2.53P ) between which merger rates are 
calculated. 

Perfect agreement between the loss and gain calculations would result in (Lg) = {Gg). As 
can been seen from the figure the fit is not perfect but is close for the most part, with the 
largest deviation being at late times in the E2A model. Comparing the three E4B runs 
shown, it is clear that including a minimal stellar bar in the model results in a substantial 
increase in the merger rates, but using an Arches-style IMF produces much larger merger 
rates from the start which also increase more over the simulation time; this will be explored 
more in Chapter [5j 

These results, and similar ones for other models, are taken to validate the overall method 
of calculating Lg and Gg. However as the fit between mass loss and gain rates is not 
perfect, at each timestep the calculated gain rate Gq{I,J) is multiplied by a correction 
factor {Lg) / (Gg) , i.e., Gg{I,J) is normalized so that the effective value of the net rate 
of mass gain due to mergers (Gg) agrees with the average rate of mass loss (Lg). {Gg is 
adjusted to agree with Lg instead of the other way around because Lg is a much simpler 
calculation and in practice shows less fluctuation over simulated time.) 

4.6 Testing the Differencing Scheme 

In order to test the differencing scheme analytic solutions of the Fokker-Planck equation 
()3.4p for both constant and varying Dij were employed. Both solutions have the form 
f{x,y,t) = j^^kxx+kyy+ujt ^ constant-coefficient case, the dispersion relation thus 

produced is w = ^Dxxk"^ + Dxykxky + ^Dyyky. For an analytic solution with non-constant 
diffusion coefficients, if one takes D^y = 0, Dxx{x) = Pxe~^'^^ and Dyy{y) = Pye"^^^ with 
Pi = const, then uj = and a static system (^ = 0) results. 

Figure 14.121 shows that the numerical boundary condition of ()3.38dp tracks the true solution 
closely, only failing to keep up on the edge at which the per-timstep change is greatest - this 
is not surprising given that ()3.38dp implements a natural spline to deal with the edge effects 
and will underestimate changes for which the second derivative is far from zero. This is not 
a large concern in practice, as tests show the discrepency is worst when the "crossterm" 
diffusion coefficient Dxy is of similar magnitude as the smaller of Dxx and Dyy (as is the 
case in Fig. I4.12p : in the actual model Dxy is typically much smaller than either Dxx or Dyy 
(although it does become of similar size when all three are close to infinitessimal) . Also, in 
the simulations the Dij don't increase exponentially near the edges, meaning the natural 
spline is a better approximation than it is in these tests. 

Figure 14.131 is similar to Fig. 14.121 except that it plots the constant-/ case and not the 
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Figure 4.11: Comparison of stellar-merger loss and gain terms Lg and Gg for the primary 
models studied, averaged over actions (/, J) and stellar-mass type. Units are arbitrary, but 
consistent within a given model. In general {Lg) and {Lg) both increase over time, so the 
progression in this plot is left-to-right for a given model. Full results for these simulations 
are given in Chapter [5l 

constant-coefficient one. 

It can been seen that use of the numerical boundary condition of (|3.38ap produced a much 
poorer match to the true solution. Use of any of the other conditions of ()3.38p resulted in 
similar or worse matches. Similarly, not using a time-splitting scheme was unsatisfactory, 
as shown by the dotted line in the figures. The Courant stability-criterion values |84] for 
the numeric tests plotted were all ^ = ^, but doubling the timestep and so doubling ^ made 
no difference in the time-split case, although it did make the non-timesplit case even worse 
than that shown. Increasing the ^ further (i.e., past ~ 1) and/or taking a much larger 
number of timesteps did start to produce larger discrepencies towards the edges where / 
is largest, presumably due to increased inaccuracies caused by the natural-spline numeric 
boundary condition more severely underestimating the second derivative of / there. 

4.7 Model Parameters 
4.7.1 Grid size: action space 

In order to have confidence in the model's results we must show that the simulations are 
stable against changes in parameters which are purely numerical: grid sizes, timestep sizes 
etc. Figure 14.141 shows a representative comparison of using different sizes of action-space 
grids. In this and several other tests, a 40x40 grid was found to give results similar to 
those of larger grids. (In most 38x38 grid was also sufficient, but anything smaller 

would deviate.) For the full simulations discussed in Chapter [5] a 40x40 grid was used unless 
otherwise noted. 
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Figure 4.12: Comparison of different finite-differencing and boundary-condition schemes. 
Tlie solid line is the analytical solution. The long- and medium-dashed lines are numerical 
solutions using ()3.38dp and (|3.38a|) respectively. The dotted line is a numerical solution also 
using (]3.38dp but in which (|3.35p is solved all at once, i.e., without the operator-splitting of 
(j3.36p . The dot-dashed line is the t = starting solution. The numeric solutions are shown 
after 30 timesteps. (Doubling At and halving the number of timesteps produced curves 
indistinguishable from those shown for the ()3.38dp and (|3.38ap cases, and a similar curve 
for the no-operating-splitting case that was only slightly different.) This figure plots a slice 
midway through the range of y used. 
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Figure 4.13: As in Fig. 14.121 but for the constant-coefficient, static-solution case (thus the 
analytic solution is identical to the t = curve at all times). Here f{x,y) is plotted for a 
slice at constant x. 
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Figure 4.14: Comparison of results of using different-size {I, J) grids. The break in the 
otherwise smooth 34x34 curve is a spurious artifact of the output-calculating procedure 
and does not affect the model's further evolution. All other parameters are the same for 
each run, and each is constrained to having the same grid as the others at each timestep. 
(These runs used grids with 50 points, a maximum / value of 9 and a 3-component 
Miller-Scalo-style IMF.) 

Choosing an even larger action-space grid is possible, but at a considerable expense of 
computer time: in addition to simply having more gridpoints on which all calculations need 
to be performed, in practice a larger grid also required smaller timesteps in order to remain 
stable. 

4.7.2 Grid size: radial coordinate 

The situation for the grid in the radial coordinate is similar to that for (/, J) in that a 
too-course grid produces evolution that is unstable to accumulated errors; however using a 
too-fine grid results in noisy output which causes difficulty for the model's mechanism of 
solving for the new potential <I>(r). In practice a 50-point grid produced a good balance 
between these two effects, although sometimes at the cost of an artificially slower evolution 
of the model as compared to finer grids. 

Figure 14.151 shows a comparison of use of 44-, 50- and 60-point grids. The 60-point case 
evolves slightly more quickly, as there are more gridpoints near the dense center of the 
cluster, but it aborts after only a few timesteps. Unlike the (/, J) grid, which is fixed, in 
the simulations described in Chapter [5] the grid is allowed to dynamically self-update 
at each timestep, as described in ^2.4.31 In practice the dynamic updating of the grid 
helps the finer-grid case adjust better than shown here: for the runs shown in Fig. 14.151 
the dynamic updating was disabled in order to allow for a direct comparison between grid 
choices. Even finer 66- and 74-point grids were also tried and gave results similar to the 
60-point case; a 99-point grid was found to require too-small timesteps with no benefit. 
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Figure 4.15: Comparison of runs with grids containing 44, 50 and 60 points; other 
parameters (timestep, action-space gridsize etc.) were identical for each of the 3 cases. 
In order to do a direct comparsion, each run was prevented from adjusting its grid 
dynamicahy. These runs were for model E4B with /max = 4 and an initial rotation parameter 
A = 0.05. The mass spectrum was 90% IM© and 10% 2Mq. 



Courser grids of 40 or 36 points could not track the center of the cluster with sufficient 
resolution. 

A somewhat more typical set of simulations is shown in Fig. 14.161 in which a broader 
range of stellar masses is used, a stellar bar is included and the grid is allowed to update 
itself dynamically at each timestep. This results in a much smoother evolution in which 
the 40, 50, 60 and 74-point grids give almost identical results, although in this case the 
simulation using a 74-point grid continued longer than the others. Note that the presence 
of the stellar bar and of a subpopulation of higher-mass stars caused much more rapid rate 
of change of the central density compared to that of the previous example. 

Given these findings a 74-point dynamically-updating grid was the default choice. In 
some simulations a 60-point grid was used; those cases are noted. 



4.7.3 Timestep size 

The choice of timestep size necessarily varied with what cluster parameters were used as 
an initial condition; a typical choice was a fraction < 0.1 of the cluster's initial central 
relaxation time trq{0) as given in ^3. It At < 0.1trq(0)- (Stellar mass ruq here is that 
of the dominant stellar mass in the cluster.) The goal was to maximize the amount of 
evolution per timestep in order to avoid inordinately long computation times, while still 
having confidence in the run's results, as demonstrated by selectively performing similar 
runs with smaller timesteps and obtaining similar results. 

Figure 14.171 shows that the simulation model is stable against different sizes of timestep. 
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Figure 4.16: Comparison of runs with grids containing 40, 50, 60 and 74 points using 
a 40x40 (/, J) grid, /max = 3 and identical timesteps in each case. Dynamic grids were 
enabled, but were similar in all 3 cases. These runs were for a 7-mass component cluster 
(mass range 1 — 32M0 in 7 bins) with a moderate amount of rotation (Aq = 0.015), Kroupa- 
style IMF and with 1% of the cluster mass in a stellar bar. 

although for different choices of timestep size the model may abort much earlier than for 
others (usually due to lack of numerical convergence when iteratively solving for the new <I> 
or p). When this happens, sometimes the run can be stably restarted from a slightly earlier 
time value but with a smaller timestep, sometimes not; for consistency all simulations 
presented are single runs starting at t = 0. 

Figure 14.181 displays a set of runs more representative of the actual simulations presented 
in Chapter m a wider range of masses is employed with a realistic IMF, and the grid is 
dynamically updated as the system evolves. Here the largest choice of At, while smooth, 
clearly does not track the progression of the central density as accurately as do the smaller 
timestep choices. The similarity of the two smaller-timestep runs is interpreted as indicating 
a reliable outcome, with the extra rise of p(0) at t = 4 — 5 Myr for the At = 0.38 Myr run 
taken to be a spurious artifact of the output-calculating procedure as seen earlier. (Note 
that in any case for the upper end of this range of masses the main-sequence lifetime is 
only 3 Myr [55].) The full span of the largest-At run is shown in Fig. 14.191 for comparison, 
and shows that the At = 0.38 run becomes possibly unstable around t ^ 8 Myr; if reliable 
results near or after t ~ 7 Myr had been required, it would have been necessary to attempt 
another smaller-At run as confirmation. 



4.7.4 Number of expansion terms 

The main physical quantity that is approximated by a series expansion is the gravitational 
potential ^'(r), which is expressed in term of spherical harmonics, as shown in (|3.14p . This 
results in the diffusion coefficients also being representated in the form of a series, as in 
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Figure 4.17: Comparison of the use of different-size timesteps. Runs are with r^-grids of 
50 points; the At = 2.4 case continued until t ~ 73, as can be seen in Fig. 14.151 Initial 
central relaxation time was trq{0) = 25.6 Myr. 
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Figure 4.18: Comparison of the use of different-size timesteps. Runs are of model E4B with 
an r^-grid of 74 points, a 40x40 (/, J) grid, a maximum / of 2 and an initial rotation given 
by Ao = 0.05 and a Kroupa IMF with 10 stellar mass bins in the range 1 — 125M0. Initial 
central relaxation time for the IMq stars that comprised 51% of the cluster was tr-i(O) = 158 
Myr; for comparison the 25M0 stars' mass fraction was 2.4% with tr7(0) = 0.23 Myr. 
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Figure 4.19: Similar to Fig. 14. 181 but showing the full range of the At = 1.52 run. 



(j3.23p ; physically, cutting off the series at a certain maximum spherical harmonic index 
I = ^max implies only summing over orbital resonances of indices equal to or lower than 
^max- Figure [1.201 shows comparisons of different choices of /max value. The cases with /max 
of 3 {i.e., up to the octopole term in the potential) and 4 do not capture the full effect of 
the dynamical friction, which can be seen as starting to converge around /max of 5. (The 
choices /max =3, 4, 5 and 6 correspond to a total number of expansion terms of 4, 6, 9 and 
12 respectively when all allowed values of index m are included.) 

When a strong stellar bar perturbation is incorporated in the potential, the situation can 
change dramatically. As comparing the timescale of Fig. I4.21l with that of Fig. 14.201 shows, 
the bar potential dominates the dynamical friction. Physically this is because the bar is a 
bulk perturbation whose total contribution to the potential is squared (c/. (13.34p ) whereas 
"field" (i.e., non-bar) stars only enter individually into the factor in (|3.23p before having 
their contributions summed over. Because the bar is by construction quadrupole-only and 
its dynamical friction timescale is so much shorter than that for field stars, higher-order 
resonances become almost irrelevent to the overall cluster evolution. So, unless otherwise 
stated all simulations which incorporate a bar perturbation have /max = 2 or 3, which results 
in only two or four terms in the expansion {e.g., the two terms m = ±2 if /max = 2). 

For the "production" simulations presented in Chapter [5] it was found that the choice of /max 
needed to be determined on a case-by-case basis. Figure 14.221 shows the behavior of one of 
the main models used, model E4B with a Kroupa IMF and a stellar bar with mass fraction 
of 1%. In this case /max = 2 gave a fairly steady evolution; higher values of /max hinted at 
a faster increase in central density p(0) but were not sufficiently numerically stable to be 
considered robust; this /max = 2 was used, understanding that it may not quite capture the 
full progress of the system. 

Contrasting this is the case shown in Fig. I4.23| of model E2A with a Kroupa IMF. Here it 
is clear that /max = 5 and /max = 6 give almost identical results, while /max < 4 misses a 
considerable degree of the increase in central density. However, the /max ^ 5 runs terminate 
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Figure 4.20: Demonstration of the use of an Zmax value of 3, 4, 5 and 6 in the diffusion 
coefficient expansion ()3.23p (and by extension in the expansion of the potential ()3.14p ): this 
corresponds to a total number of expansion terms of 4, 6, 9 and 12 respectively. This plot 
is of a rotating 2-component cluster similar to those of Figures [4.151 and 14. 171 Each run was 
constrained to use the same timestep and the same r^-grid as the others, to allow for direct 
comparison. 
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Figure 4.21: Comparison of different numbers of expansion terms and gridpoints for 
the case including a bar perturbation. The 50-point, quadrupole-only run continued on a 
fairly linear path (not shown) until it reached p ~ 2 x 10® at t = 0.2. This run is for a 
2-component cluster with 10% 2Mq and 90% IMq stars, 4% of which comprised the bar. 
The rotational parameter was A = 0.075. A run with l^^x = 4 (i.e., 6 expansion terms) was 
nearly identical to the one shown using Z^ax = 3, although it ended earlier. 
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Figure 4.22: Comparison of different numbers of expansion terms for a "production" run: 
model E4B with a Kroupa IMF, initial rotation parameter Aq = 0.05 and a stellar bar mass 
fraction of 1%. The stellar mass range is 1 — 125 Mq. 
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Figure 4.23: Comparison of different numbers of expansion terms for a "production" run: 
model E2A with a Kroupa IMF, initial rotation parameter Aq = 0.05 and a stellar bar mass 
fraction of 1%. The stellar mass range is 1 — 125 Mq. 
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Figure 4.24: Comparison of central density for runs in which the stellar mass spectrum is 
split into the number of geometrically-spaced mass bins shown. All are model E4B with 
50 points in the grid, a 34x34 (/, J) grid, and an overall mass range of 1 — 32Mq. The 
IMF is Salpeter-like, but with the uppermost mass bin artifically overpopulated in order to 
accentuate any differences in the cases. 

quite a bit earlier than the lower-Z^ax simulations, and require considerably more compu- 
tational time (weeks compared to days) to even get as far as they do. (Taking /max — 4 
produces a total of 6 terms in the expansion, while /max = 6 gives 12 terms). So in the 
case of model E2A/Kroupa an /max of 4 was used, but only to give an indication of how the 
system evolves and not to show precise quantitative results. 

The above two examples are the extreme cases of those studied; in general the bar perturba- 
tion proved dominant enough that /max = 2 or 3 sufficed to describe the system's evolution 
while still allowing for numerical stability and adequately-short computation times. 

4.7.5 Mass Spectrum: Discretizing the Initial Mass Function 

The initial mass functions described in ^2.6.31 are continuous functions of the individual 
stellar mass m; in order to separate a given IMF into discrete bins of discrete stellar mass ruq 
the IMF is simply integrated over the range of masses being considered, with the boundary 
between bins being the midpoint between them. The lowest mass bin's lower limit is taken 
to be ^mi, and the highest-mass mg bin is given the same width as it would have if there 
was an mg+i bin with the same scaling. As the spectrum of mass bins is geometrically 
increasing by a factor ruq^i < 2mq, this prescription artifically biases the numerical IMF 
slightly towards lower masses as compared to the analytical form on which it is based. 

Figure 14.241 shows the evolution of the central density for test runs of model E4B with a 
stellar mass range between 1 and 32M0 and a Salpeter-like IMF in which the 32M0 bin 
has been overpopulated, so that high-mass stars drive the system's evolution more strongly 
than would occur naturally. Figure 13 . 2 5 1 simlarlv shows results for the same runs, but for the 
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Figure 4.25: Similar to Fig. 14.241 but showing the central density of the uppermost (mg = 
32M0) stellar mass bin only. 
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Figure 4.26: Comparison of central density for runs in which the stellar mass spectrum is 
split into the number of geometrically-spaced mass bins shown. All are model E4B with 50 
points in the grid, a 34x34 (/, J) grid, and an overall mass range of 1 — 32M0. The IMF 
is Kroupa-like, similar to that used in the later full simulations. Rotation and a stellar bar 
are included, with a bar mass fraction of 1% and an initial rotation parameter of Aq = 0.016. 
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Figure 4.27: Similar to Fig. 14.261 but with stellar mergers included in the simulation. (The 
9-bin case shows shghtly lower central density than the others due to it following a somewhat 
different series of values for the innermost point of its dynamic grid; this was unusual 
and did not occur during the full simulations presented in Chapter [5] except where noted.) 

32Mq stellar mass bin only. Although the curves of the central density in this artifical case 
are somewhat rough, a choice of 7 or 8 mass bins seems to best show the full evolution of 
the system while avoiding the numerical instability displayed by the 9-bin run. For a range 
of 1 — 32Mq, 7 bins corresponds to a mass ratio between adjacent bins of 1.78, and 8 bins 
gives a ratio of 1.64. The full simulation runs shown in Chapter [5] used adjacent-bin mass 
ratios of 1.71 for the Kroupa IMF and 1.675 for the Arches-style IMF. This is similar to and 
consistent with the findings of Amaro-Seoane |4l] who found that an average mass-bin ratio 
of ~ 1.72 was sufficient to model the Fokker-Planck evolution of clusters of 0.2 — IOOMq 
stars which contained already- formed or primordial massive central objects. 

Figures [4. 261 and 14. 271 show the central density evolution using a more-realistic Kroupa-style 
IMF and incorporating a stellar bar and a moderate amount of rotation. It can be seen 
that both with and without stellar mergers, the choice of number of stellar mass bins does 
not strongly affect the simulation results; as seen in Fig. 14.271 the dynamic grid can 
sometimes have a greater, but still small, effect. 

4.8 Binary Heating 

As binary stars' orbits harden as a result of encounters with other stars, the system is 
effectively heated by the energy tranferred from the binary system to the cluster as a whole, 
and potentially disrupted somewhat. The physics of binary heating is not incorporated in 
the simulation code; however, an upper limit on its possible effect was calculated for a 
selection of models using a direct application of the binary-heating formulation of Quinlan 
and Shapiro [7j. The extreme assumption that all stars in the cluster are in binaries and 
that all binaries harden forever (i.e., none are soft binaries that do not harden, and none are 
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Figure 4.28: Ratio of maximum possible amount of binary heating Kb to overall (gravi- 
tational + kinetic) system energy \Etot\ for a variety of models. All used the Arches-style 
IMF except for model "E2A/Kroupa", and all had an initial Plummer-sphere distribution 
except for "G2A", which started as a 7 = sphere. A logarithmic scale is used so that 
different models' timescales can be plotted together; exactly flat lines at the start of each 
model's plot indicate the first timestep in the calculation, before which binary heating is 
zero. 



themselves disrupted by interactions with other stars) is used to give an absolute limit on 
how much binary heating could possibly have occured, had it been included in the dynamical 
calculations. 

The results are shown in Fig. 14.281 from which it can be seen that at all times (except 
at the very end of one simulation) the upper limit of possible binary heating remains a 
small fraction, 10^^ or less, of the overall system energy. Given that not all stars can be in 
binaries, that not all binaries will harden and heat the system, and that some binaries will 
be disrupted by encouters and by internal merging, the assumption that binary heating can 
be neglected holds. 
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Chapter 5 

Results 



5.1 General Considerations 

5.1.1 Initial Models 

The progression of initial potential-density models studied starts with the E2A and E2B 
Plummer spheres which, as listed in Table 12.11 and described in ^2.6.H can represent the 
nuclei of dwarf elliptical or bulgeless spiral galaxies, as well as very dense globular clusters. 
Model E4B is somewhat more massive and more dense with a higher initial velocity disper- 
sion, along the lines of the core of a giant elliptical galaxy. Model E4A starts out even more 
dense than E4B, but was found to be numerically unstable and so could not be simulated. 

An alternative to the Plummer sphere is an initial "7 = sphere", which can represent 
galactic spheroids. When mapping the masses and core radii of the above Plummer models 
onto 7 = spheres, only the resulting model G2A (analogous to Plummer model E2A) 
has an initial central density and velocity dispersion within the ranges required for mod- 
eling dense astrophysical systems without assuming an unreasonably high central density. 
However, a model G3C, which starts with the highest initial density studied, was run for 
comparison purposes, as was intermediate model G3A; these are described later in this 
Chapter. 

For each of the various potential-density models the distribution of stellar masses is given 
by the initial mass function (IMF). As developed in ^2.6.31 the base IMF used here is the 
Kroupa IMF, which is a good fit for the observed galactic stellar population in general. The 
second is an Arches-style IMF which has recently been determined to specifically describe 
the stellar distributions in dense clusters observed at the centers of galaxies, i.e., the objects 
being simulated here. Unsurpringly the Kroupa IMF is more weighted towards lower-mass 
stars while the Arches IMF has a fiatter distribution of masses, but still with an upper 
cutoff of m < 150Mq. 

5.1.2 Physical Effects 

For sufficient buildup of the cluster's density to occur so that a massive object can form 
as the product of collisional mergers of stars, rotational support against the infall of stellar 
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material to the central region must be overcome. As described in Chapter [H given infinite 
time such collapse would be expected even in a nonrotating system if it is sufficiently dense. 
However, if it is to occur before the majority of stars reach the end of their main-sequence 
lifetimes - only 3 Myr for the largest stars in the IMFs considered - something in addition 
to interactions of individual stars is expected to be required. Fortunately, even though 
they were not intentionally constructed to be so, when given a reasonable value for the 
rotation parameter A the models chosen all satisfy the criterion for being unstable against 
the creation of a bar-like perturbation in the stellar potential; this stellar bar can (hopefully) 
then transport angular momentum in bulk from the system's center to its outer regions. 

Thus the problem is essentially a race against time: can gravitational effects, aided by 
the transport of angular momentum afforded by a stellar bar, allow the system's central 
regions to condense quickly enough that collisional mergers can in turn produce an object of 
M > 25OM0 which will evolve into a massive seed black hole, all before the stars outlive their 
main-sequence lifetimes. Simulations for each model and IMF are performed incorporating 
all of: an overall rotation of the system, a stellar bar, and collisional mergers of stars. In 
order to examine the various effects some simulations are also presented in which the system 
is not rotating, it does not contain a stellar bar, and/or the stars are not allowed to collide 
and merge to produce larger stars. 

After a brief summary of what constituted successful simulations and which models yielded 
them, the remainder of this chapter reports results for each choice of initial model and IMF. 

5.1.3 Overview 

Of the five possible initial models whose density-potential pairs are listed in Table 12.11 and 
described in ^2.6.11 (or eight, when models G2A, G3A and G3C from ^2.6.11 are included), 
only some yielded simulations which satisfied all the criteria for stability. The criteria 
included: 

• giving consistent results with different choices of timestep size At; 

• there existing a value of the cutoff /max of the potential expansion high enough that 
terms / > /max did not contribute significantly (as indicated by the central density's 
increase with time), but low enough to avoid the numerical instability higher values 
of /max exhibited; and 

• giving consistent results across the Sun, Alpha, and AMD/x86 computer architectures. 

The models which met the above criteria are listed in Table 15.11 For each model a possible 
astronomical system to which it best corresponds (if any) is listed; this expands on the 
descriptions given in §2.6.1} which did not yet take choice of IMF into consideration. Table 
15. il can serve as a convenient reference for the model- specific sections that comprise the 
remainder of this chapter. 

Simulations with a Kroupa IMF used a mass range of 1 — I25M0; the Arches IMF has a 
smaller contribution from low-mass stars (c/. ^2.6.3p and so employed a range of 2 — 125Mq. 
The set of possible variations of each model include: whether the system has initial rotation 
(with the rotation parameter set near the canonical value A ~ 0.05 as described in ^2.6.2p 
or not (A = 0); if rotating, whether there is a stellar bar or not; and whether stellar mergers 
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Initial Profile Model Set IMF Corresponding Astronomical System 





E2A/E2B 


Kroupa 


dwarf elliptical or bulgeless spiral galaxy nucleus 


Plummer 




Arches 


nuclear cluster 


sphere 


E4B 


Kroupa 


core of a giant elliptical galaxy 






Arches 


none - inappropriate IMF 




G2A 


Kroupa 


galactic spheroid or spiral bulge 


7 = 




Arches 


none ~ inappropriate IMF 


sphere 


G3A 


Kroupa 


galactic spheroid or spiral bulge 




G3C 


Kroupa 


none - initial central density too high 



Table 5.1: Models which yielded stable results, along with possible analogous astronomical 
systems. 



are allowed or not. For both the Kroupa and the Arches-style IMF, two models produced 
stable results in at least some cases: E2A and E4B. The E2B model was stable as well but 
only when using the Arches IMF. The results for each model are described in detail below, 
and are summarized in Tables 15.21 and 15.31 for the "full" runs which included rotation and 
stellar mergers. 

5.2 Plummer-sphere Models 

5.2.1 Model E4B, with Kroupa IMF: Core of a Giant Elliptical Galaxy 

With Mtot = 3.1 X WMq and an initial p{0) = lO^Mo/pc^, Model E4B was the largest 
overall system studied with the highest central density. (Model E4A had an initial p(0) = 
3 X IO^Mq/pc^ but did not yield stable simulations.) Figure [5TT] shows the change in central 
density p{0) with time for the four cases: (1) nonrotating, (2) rotating, (3) rotating, with 
a stellar bar and (4) rotating, with a stellar bar and star-star mergers. The stellar bar 
clearly influences the system evolution greatly as it produces an increase in p{0) of more 
than an order of magnitude in the 3 Myr main-sequence lifetime of the largest stars in the 
system, compared to an increase of < 5% for the cases without a stellar bar. Contrasting 
this is the effect of stellar mergers, which is minimal: at the end of the simulation that 
included mergers, < 0.1% of the lowest-mass stars had undergone a collision and merger, 
while fewer than 30 of the highest-mass (I25M0) stars had been created from collisional 
mergers - although this increase is still visible in Figure [5^21 Table [52] shows that G250, the 
calculated rate of producing 250Mq stars by stellar mergers, is effectively zero even after 3 
Myr; similarly, the production rate for all stellar masses larger than 2OOM0 was also found 
to be zero. 

Looking at the system as a whole, Fig. 15. 31 shows the range of p{r) over the full radial extent 
of the stellar cluster at the starting time, the midpoint, and the end of the simulation. As 
the core contracts and gets more dense, the outer regions lose stellar density accordingly. 
Interestingly, when the model is run without including a stellar bar no such change in 
density is observed in the outer regions; this behavior is seen in general for all the Plummer- 
sphere models studied and is attributable to the stellar bar transporting angular momentum 
outwards, which rarifies the outer regions while allowing the inner regions to lose rotational 
support and contract. 
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Model /max bar p(0) at to p(0) at h h ^250(^0) ^250(^1) ^(^1) S'ih) 



E2A 


4 


5.3% 


2.5 X 10'' 


3.0 X 10« 


1.6 


0. 


0. 


0.82 


0.83 


E4B 


2 


none 


6.3 X 10'^ 


7.0 X 10'^ 


2.5 


0. 


0. 


1.09 


1.11 






none 




7.2 X lO'^ 


3.0 


0. 


0. 


1.09 


1.11 






10% 


;5 


3.8 X 10^ 


2.5 


0. 


6 X 10-5 


1.006 


1.02 


G2A 


4 


none 


1.7 X 10^ 


1.3 X 10« 


.42 


0. 


0. 


0.95 


0.96 


G3A 


3 


7.7% 


4.4 X 10^ 


1.1 X 10^ 


.46 


0. 


0. 


0.87 


0.87 


G3C 


2 


7.7% 


1.4 X 10« 


3.0 X 10'^ 


.16 


0. 


0. 


1.04 


1.04 



Table 5.2: List of models which yielded stable simulations using the Kroupa IMF, with 
central-density and merger-rate results shown at various times ti in each simulation (usually 
the endpoint, or chosen for comparison with another). Times are in Myr, central densities 
in Mq/pc^, and each simulation had a rotation parameter A ~ 0.05. 6*250 is the calculated 
rate at which 250Mq stars would be expected to be produced via collisional mergers of 
125Mq stars, in Myr~^. S is the factor by which the central density of the most-massive 
stars increased as a result of mass segregation, relative to the overall increase in central 
density; S' is the same factor but also including the effects of stellar mergers. For the G2A 
and G3C models, the presence or absence of a stellar bar had little effect on the system's 
evolution. (Model G3A was not run without a stellar bar but it was also expected to have 
little effect.) Note that the central densities listed here are for the initial models with 
rotation, and so are lower than those given in Table [2Tl Starting time to = 0. 



Model 


/max 


bar 


/)(0) at to 


p{0) at h 


h 


6*250 (^o) 


G*25o(il) 


S{h) 


S'ih) 


E2A 


3 


5.3% 


2.5 X 10'^ 


3.7 X 10« 


.63 


3.9 


6.5 


1. 


1.004 


E2B 


2 


none 


8.2 X 10" 


1.4 X 10^ 


3.0 


2.2 


2.7 


1. 


1.01 






5.3% 


1? 


6.1 X 10^ 


3.0 


2.6 


3.3 


1. 


1.01 




3 


17 


11 


9.4 X 10'^ 


3.0 


2.6 


3.7 


1. 


1.01 






8% 


1? 


5.3 X 10*^ 


2.7 


2.3 


8.3 


1. 


1.006 


E4B 


2 


8% 


6.3 X 10*" 


1.0 X 10*^ 


2.2 


35. 


73. 


1. 


1.015 




11 


10% 


1; 


1.5 X 10^ 


2.0 


33. 


69. 


1. 


1.025 


G2A 


3 


none 


1.7 X 10 '' 


1.3 X 10« 


.07 


1.2 


(1.9)* 


1. 


1. 




51 


7.7% 


1.7 X lO'^ 


1.3 X 10® 


.07 


1.1 


1.6 


1. 


1. 



Table 5.3: Similar to Table [5^21 but for simulations using an Arches-style IMF. [*] indicates 
that the value of ^250(^1) calculated for model G2A without a bar is considered to be 
numerically suspect. 



Model 


/max 


bar p{0) at to 


^(0) at ti 


h 


G*25o(io) ^250(^1) 


s{h) 


S'ih) 


E2B 


2 


5.3% 8.2 X 10*^ 


2.9 X 10® 


4.3 


2.6 5.1 


1. 


1.01 




3 


71 17 


4.9 X 10® 


4.3 


2.6 8.7 


1. 


1.01 



Table 5.4: Addendum to Table [5T3l showing results of simulations allowed to proceed beyond 
the 3 Myr lifetime of the largest main-sequence stars. 
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Figure 5.1: Central density v. time for model E4B using the Kroupa IMF with a stellar 



mass range of 1 — 125M0 and lr_ 



2. Note that the nonrotating (Aq = 0) case has an 



initial p{0) slightly higher than the cases that include rotation; this is a product of the 
method for introducing rotation into the initial distribution function of the system. 
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Figure 5.2: Similar to Fig. 15.11 but showing only the highest-mass stars (mg = I25M0). 
The upturn at the end of the "with mergers" curve is for the final timestep only and is a 
result of a numerical instability which casued the simulation to end early. 
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Figure 5.3: Density p{r) v. radial distance r from the cluster center for the start, midpoint 
and end of the simulation for model E4B with a Kroupa IMF. Stellar mergers were enabled, 
as was a 10% stellar bar. Rotation parameter A = 0.05 and /max = 2. (This corresponds to 
the uppermost plot of Fig. 15. li ) Note that the decrease in p{r) with time for large r is only 
exhibited when a stellar bar is present. 



The quantity S in Table 15.21 gives the relative increase in the central density of the highest- 
mass (125 Mq) stars relative to that of the overall system, and so indicates how much 
mass segregation has occured. S' is a similar measure but also includes the effects of stellar 
mergers on the increased density of high-mass stars. For model E4B without a stellar bar 
a moderate (< 10%) relative increase in high-mass stars is seen, with mergers contributing 
another couple percent. (The nonrotating A = case, plotted in Figures 15.11 and 15.21 but 
not listed in Table [5^21 displayed similar behavior.) In this measure the presense of a stellar 
bar appears to strongly damp the relative mass segregation to a large extent, but that is 
merely due to how much more the overall central density increases with the bar (by a factor 
of ~ 6 with the bar as compared to ~ 1.1 without). Thus a small relative mass segregation 
in the with-bar case still corresponds to a larger absolute increase in the central density of 
the highest-mass stars. Still, it also means that the bar is dominating the dynamics and 
the different stellar-mass populations evolve less distinctly than in the no-bar case. 



5.2.2 Model E4B, Arches-style IMF 

When studying model E4B with the Kroupa IMF replaced by an Arches-style IMF the 
situation changes quite a bit. As shown by Figures 15.41 and 15.51 both the overall central 
density and that of the largest-mass stars increase much more rapidly even in simulations 
without a stellar bar: by 3 Myr, p{Q) is up by a factor of > 3 — 4 for the nonrotating and 
rotating cases, even without the presence of a stellar bar. More dramatically, the systems 
which included a bar perturbation showed an increase in p{<d) of ~ 20 x or more, and in 
a shorter amount of simulated time (by 2 Myr when using a bar of strength 10%, by 2.2 
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Myr for an 8% bailll), after which the simulations were not able to track the evolution 
stably. This even more rapid increase is consistent with the Arches IMF being weighted 
more towards larger stellar masses than is the Kroupa IMF. The figures show that unlike in 
the Kroupa-IMF case, with an Arches IMF the stellar system achieves core collapse before 
the 3 Myr stellar-lifetime limit, with the central density p(0) increasing exponentially: the 
run with a 10% bar and mergers was able to track somewhat into the core collapse stage 
before aborting, while the other with-bar runs reached the turning point of p{0). Once core 
collapse begins, the central core of the cluster decouples dynamically from the outer regions 
- a situation the Fokker-Planck model is not constructed to deal with. 

The stellar merger situation is also qualitatively different with an Arches IMF: as seen in 
Table [5^ even at t = the rate at which 25OM0 stars were being produced from collisions 
of 125M0 stars was ~ 33Myr~^, increasing to ~ 69Myr~^ by the end of the simulation 
with a 10% stellar bar. Thus the Arches IMF was sufficient to produce 250Mq stars from 
the outset, although the presence of a stellar bar enhanced the production. It was also 
found that by t = 2.0 Myr, approximately 500 stars of mass 125M0 had been produced by 
collisional mergers of lower-mass stars. 

Interestingly, the flatter Arches IMF seems to strongly damp any mass segregation when 
compared to the Kroupa-IMF case above: the relative increase in central density Pq(0) of 
the highest-mass stars tracks that of the overall system p{0) to several significant figures, 
regardless of the presence or strength of the stellar bar. And again, the nonrotating A = 
case plotted in Figures 15.41 and 15.51 (but not listed in Table 15. 2p displayed a similar complete 
lack of mass segregation; 5 = 1 for all E4B/ Arches cases. There is only a slight relative 
increase in pq(0) due to stellar mergers, as shown by S'{ti) = 1.025. 

Figure ISlBl plots the stellar mass density p(r) from the center of the cluster outwards. As in 
the Kroupa-IMF case, the stellar bar transports angular momentum from the inner regions 
to the outer, allowing core contraction and a corresponding decrease in stellar density at 
large r. A similar simulation but which does not include a stellar bar is plotted in Fig. 15.71 
Without a bar to transport angular momentum outwards, there is very little decrease in 
stellar density in the outer parts of the cluster. Note that Figures 15.61 and 15.71 are plotted 
on the same scale for easier comparison. 

Other quantities: bar composition, velocity dispersion and relaxation time 

For the "with bar" runs presented in this chapter, the default was to populate the bar with 
stars from the lowest-mass stars considered in the simulation, e.g., 2Mq stars when using 
the Arches IMF. As a test, the E4B/Arches "with bar" simulation was also run with a bar 
comprised of IMq stars. It gave identical results to the default case, which is expected 
since for purposes of its perturbing potential the bar is treated as a single bulk object and 
not collection of individual stars. 

In addition to stellar density, several other quantities are tracked by the simulation code. 
Two of the most interest physically are the velocity dispersion ai and the relaxation time 
tr- Model E4B/Arches can be used as an example to show the results for these properties, 
which evolved similarly for all models. 

^In model E4B the strength-10% bar case used the default 1% mass-fraction bar, as described in i]3.6l the 
strength-8% bar corresponded to a mass fracton of 0.8%, below which the code had difficulty converging. 
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Figure 5.4: Central density v. time for model E4B using the Arches IMF with a stellar 
mass range of 2 — 125M0 and /max = 2. 
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Figure 5.5: Similar to Fig. 15.41 but showing only the highest-mass stars [rnq = I25M0). 
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Figure 5.6: Density p{r) v. radial distance r from the cluster center for the start, midpoint 
and end of the simulation for model E4B with an Arches-style IMF. Stellar mergers were 
enabled, cis Wets 3> 10% stellar bar; rotation parameter A — 0.049 and. /max — 2. (This 
corresponds to the leftmost plot of Fig. 15. 
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Figure 5.7: Density p{r) v. radial distance r from the cluster center for the start and end 
of the simulation for model E4B/ Arches, as well as at two similar times as those shown in 
Fig. 15.61 for comparison. Neither stellar mergers nor a stellar bar were enabled (in contrast 
with Fig. 15. 6p . Rotation parameter A = 0.049 and /max = 2. 
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Figure 5.8: Central velocity dispersion o"i(0) v. time of the lowest-mass {2Mq) stars in 
model E4B with an Arches IMF, for nonrotating, rotating and with-bar cases. (The case 
shown with a stellar bar also had stellar mergers enabled, but a run without mergers was 
similar.) 




I I I I I I I I 

0.5 1 1.5 2 2.5 3 

t [Myr] 

Figure 5.9: Median central relaxation time tr(0) v. time of the lowest-mass (2M0) stars in 
model E4B with an Arches IMF, for nonrotating, rotating and with-bar cases. (The case 
shown with a stellar bar also had stellar mergers enabled, but a run without mergers was 
similar.) 
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Figure 15.81 shows the change with time of ai (0) in the various cases studied. The systems 
which are initially rotating start at a lower value of cri(O) because the method of incorpo- 
rating rotation into the initial distribution function /(/, J, 0) necessarily "cools" the system 
by increasing the amount of order present in the distribution of angular momentum J. The 
two runs which do not include a stellar bar have a slowly increasing (Ti(0) as the system 
contracts and the central core gains entropy. Contrasting this is the case with a stellar bar, 
which collapses much more quickly and has a corresponding rapid increase in the central 
velocity dispersion. Of note, the expected initial value of o"o(0) for model E4B without 
rotation as given by Table [2m is 400 km/s; all models had a calculated initial cti(0) which 
was similarly lower than would be expected if isotropy of the velocity dispersion (as given 
by ()2.39p with (5 = 0) strictly held. To make the calculated initial cti(0) values agree with 
those of o"o(0) in Table \27i\ an anisotropy parameter 6 ~ —0.3 is required; fortunately no 
derived quantities in the calculation depend strongly on an assumption of isotropy. 

The median central relaxation time tr{0) for the same set of cases is shown in Fig. 15.91 The 
initial value of tr(0) is also slightly different than that listed in Table [2Tt which is expected 
as the full calculation of tr applies to a given test star, not to a specific point in space 
(c/. the footnote in §3.ip . Thus only a median value of tr across all stars which traverse 
the system center can be determined, and any slight error in cti(0) does propogate to the 
deduced tr{0), via its dependence on the cube of ai. This being the case, it can be seen 
that the presence of a stellar bar does cause a rapid decrease in the central relaxation time 
of the system, in this case by close to an order of magnitude within 2 Myr. 

5.2.3 Model E2A, Arches-style IMF: Nuclear Cluster 

Model E2A is somewhat different than E4B: the total mass, initial central density and 
velocity dispersion are all lower for E2A, as is the initial relaxation time (c/. Table [2T]) . 
When using the Arches IMF the simulations which did not include a stellar bar were not 
numerically consistent for different choices of the timestep size, and so only the "full" case 
that included a stellar bar and collisional mergers is shown in Fig. 15.101 

Despite these differences, the increase in p{0) by a factor of ~ 14 before the simulation ended 
was similar to that of model E4B prior to its core collapse. Model E2A did not reach core 
collapse before aborting - although for E2A the simulation ended much earlier, at t = 0.63 
Myr, consistent with the shorter initial value for t.,. and smaller ai. And similarly to as 
was found for model E4B, there was sufficient merging of 125Mq stars initially to create 
250Mq stars - in this case, at a rate of 3.9Myr~^, which increased to 6.5Myr~^ by the 
simulation's end. Approximately 15 stars of mass 125M0 had been created via collisional 
mergers of lower-mass stars as well by the simulation's end, a small number attributable to 
the short amount of simulation time. The mass segregation behavior, with no segregation 
observed {i.e., S = 1) due to dynamical effects and only a small amount (S' > 1) from 
stellar mergers, was also similar to that of the E4B/Arches model, as was to be expected 
given the short simulated time of evolution of the system. 

5.2.4 Model E2B, Arches-style IMF: Nuclear Cluster 

Model E2B, when given an Arches-style IMF, shares many properties with model E4B/Arches 
described in §5.2.21 Figure [5TT] shows that the nonrotating case evolves very slowly, as does 
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Fi gure 5.10: Central density v. time for model E2A with Arches-style IMF, /max — 

3 and a 

stellar mass range of 2 — 125M0. 
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Figure 5.11: Central density v. time for model E2B with /max = 2 and using the Arches 
IMF with a stellar mass range of 2-125Mq. Note that the main-sequence lifetime of the 
most massive stars is 3 Myr and so the region of the graph beyond f = 3 is nonphysical. 
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Figure 5.12: Detail of the first 3 Myr of Fig. KTl\ 
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Figure 5.13: Similar to Fig. 15.111 but including Zmax = 3 and an 8% bar case. All have 
A = 0.049 and do not include stellar mergers. (Simulations performed with mergers enabled 
yielded almost-identical runs of p(0) to those shown in this figure.) Note that the "5% bar" 
plot here is the same as that shown in Fig. 15.111 The gap in one plot is due to a glitch in 
the output routine, as described in Chapter [H 
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Figure 5.14: Density p{r) v. radial distance r from the cluster center for the start, midpoint 
and end of the simulation for model E2B with an Arches-style IMF. Stellar mergers were 
enabled, as was an 8% stellar bar; rotation parameter A = 0.049, and /max = 3. This 
corresponds to the leftmost plot of Fig. 15.131 The decrease in p{r) with time for large r is 
only seen when a stellar bar is present. 



the rotating case which does not have a stellar bar; in contrast, a modest 5.3% stellar bar 
induces core collapse and a ~ 60x increase in central density p{0) ~ but only at t ~ 4 Myr 
of simulation time, well beyond the 3 Myr main sequence lifetime of the largest stars. As 
seen in Fig. 15.121 the evolution prior to t = 3 Myr is much more modest. Stellar mergers 
have little affect on p{0), and Table shows that with a 5.3% bar, by t = 3 Myr the rate 
of massive star production increases modestly from G250 = 2.6 Myr~^ to 3.7 Myr~^. [i.e., 
a small number of massive stars are produced collisionally at the start, and a small number 
are produced at the end of the stellar lifetimes.) Table 15.41 shows that if the system was 
able to continue on to t ~ 4 Myr, core collapse would result in a somewhat larger increase 
in the rate of massive-star production. 

A rise over time in the rate of producting 25OM0 stars via mergers can be attributed either 
to an increase in overall density p, or to the presence of a greater number of massive stars 
due to previous collisional mergers. A comparison to determine how much each of the above 
factors contributed was performed by running a simulation in which stellar mergers were 
shut off at first and then turned on at t = 3 Myr. The result was a value of G250 — 3.4 
Myr~^ at t = 3, compared to 2.6 Myr~^ at t = 0, and to 3.3 Myr~^ at t = 3 when mergers 
were allowed from the start. Hence, at least for model E2B with a 5.3% bar, there is little 
cumulative effect of collisional mergers on the later merger rate of larger stars. 

Allowing for a somewhat stronger bar gives a significant effect: Figure 15.131 shows that 
with an 8% bar (still well within the typical range of 5 — 10%) core collapse is accelerated 
and shifted earlier to t ~ 2.5 Myr; likewise the rate of 25OM0 star production increases 
somewhat, now to 8.2 Myr~^ - still not as high as for model E4B/ Arches but notably more 
than at t = in this model. As well, ultimately ~ 60 stars of 125Mq were created by mergers 
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Table 5.5: Effect of artificially increasing the collisional merger rate by a factor of 5 in 
order to account for the effect of tidal-capture binaries, as described in ^5.2.61 Listed are 
results for model E4B with /max = 2 and a 10% stellar bar. The factor = 1 entry has slight 
differences from the values shown in Table 15.31 due to different timestep sizes being used for 
the simulations listed here. 

of lower-mass stars by the time the simulation ended at t = 2.7 Myr. And also similarly to 
model E4B/ Arches, no mass segregation is seen before the simulations terminate - which 
again was also the case for a simulation performed with A = {i.e., without rotation). 

Finally, the full range of stellar mass density p{r) plotted at various times in Fig. 15. 141 again 
shows the core of the cluster contracting and becoming denser, while the outer regions lose 
density over time. As was also seen for model E4B, a simulation performed without a stellar 
bar did not display any similar effects of angular momentum transport to the very outer 
regions of the cluster as did one which included a bar. (The non-bar case is not plotted 
here, but compares similarly to Fig. 15.141 as do Figures l5. 131 and 15.61 for E4B.) 

5.2.5 Model E2A, Kroupa: Bulgeless Spiral or Dwarf Elliptical Nucleus 

The simulations of model E2A with the Kroupa IMF behaved somewhat strangely: for very 
similar cases - e.g., two runs both with rotation and a stellar bar, but only one with stellar 
mergers (which had a small effect when using the Kroupa IMF) - the dynamic r"^ grid 
would evolve rather differently for the various cases, making direct comparisons difficult. 
Also, finding a satifactory value for the potential -expansion parameter /max 

proved elusive: 

^max = 5 appeared to be required to capture the bulk of the interaction strength, but was 
even more unstable than /max < 4. Still, as shown in Table 15.21 this model gave results 
consistent with model E4B/Kroupa: the stellar bar produced an increase in yo(0) of an 
order of magnitude in a short time compared to the main-sequence lifetime of the most 
massive stars, but the rate of collisional mergers of those massive stars was negligible. The 
"reverse mass segregation" indicated by S < 1 is unexpected. However, test runs using 
a lower value of /max = 2, while not capturing the full extent of the model's evolution, 
also showed a similar trend towards S < 1 through t = 1.6 Myr - but at later times S 
turned around and began increasing again. Even so, why there would be a delay in mass 
segregation being exhibited for this model remains unexplained. 
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Figure 5.15: Effect of artificially increasing the collisional merger rate. Plotted is model 
E4B using a Kroupa IMF with /max = 2 and a 10% stellar bar. Also see Table [531 

5.2.6 Effect of Collisional Merger Rates 

The development of the collisional merger rates in ^3.81 assumed only direct 2-body colli- 
sions. However, indirect collisions - due to the creation of tight binaries formed via tidal 
capture which then merge on a short timescale - can enhance the merger rate by a factor 
of 3 to 5 in some situations [88]. (Other factors which may enhance the expected rate 
include that low mass [5 — IOMq] stars do not to contract to their long-term main sequence 
radius immediately and so have a larger cross section initially, and a second-order effect in 
that when those lower-mass stars collide with more massive stars they form a disk, which 
increases the new star's effective cross section.) To test how an enhanced merger rate would 
affect the overall results, model E4B was re-run with a the loss and gain coefficients of 
(j3.44p and (|3.50p increased by a factor of 5, the effect of which is shown in Fig. 15.151 and 
Table [531 (Tabulated but not plotted are the results for model E4B with the Arches IMF; 
while stable, the grid of that case's factor = 5 run sufficiently diverged from that of the 
base factor = 1 simulation to preclude a direct visual comparison.) 

While in both the Kroupa- and Arches-IMF cases an enhanced merger rate resulted in a 
faster evolution of the overall system, little qualitative difference was observed. The central 
density p(0) and merger-induced mass segregation ratio S' both increased more quickly with 
the additional factor of 5 in the merger rate ~ as did the rate of production of massive stars 
^^250(^1) when using the Arches IMF. But still no massive stars were created when using 
the Kroup IMF, and core collapse did not occur any more quickly with the Arches IMF. 
As model E4B has the highest initial p(0) and is already the fastest-evolving of the various 
models studied, other models are expected to also show little qualitative change when given 
an enhanced merger rate. 
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Figure 5.16: Central velocity dispersion v. time for the lowest-mass (IMq) stars in the 
7 = sphere models with a Kroupa IMF. All have A ~ 0.05. Plots of runs with and 
without stellar mergers were identical to each other; differences within each model are due 
to the presence or absence of a stellar bar. 

5.3 7 = Sphere Models 

As described in ^2.6. H of the 7 = models which could be made using same values of 
total mass M and core radius rcorc to the Plummer-sphere models described above, only 
the "G2A" model has reasonable values for the initial central density p(0) and velocity 
dispersion cti(0). An extreme model "G3C" in which the half-mass radius was the same as 
that for model E2B was also run and is described here, as is an intermediate model "G3A" . 

Figures 15.161 and 15.171 show C7i(0) and the central relaxation time tr(0) for various cases of 
models G2A, G3A and G3C, all with a Kroupa IMF. Unlike as was seen for the Plummer 
sphere models, the stellar bar had little effect on the 7 = models' evolution, as can be 
seen from the contrast of Fig. 15.161 with Fig. 15. 9i This trend was common to all the 7 = 
sphere results and can at least partially be attributed to the short amount of simulation 
time before the models ended, as indicated by the small values on the plots' horizontal axes. 
Relaxation times were much shorter than for the Plummer sphere models, and so the bar 
had little time to transport angular momentum ~ and mergers had little time to produce 
larger stars as well. Even so, in some cases the short fr-(O) values did allow for rapid overall 
evolution as is described in the following sections for each specific model. 

5.3.1 Model G2A, Kroupa IMF: Galactic Spheroid / Spiral Bulge 

Model G2A with a Kroupa IMF showed only an almost linear increase in /3(0) before ending 
due to numerical issues; Figure [5 . 181 shows that this was the situation for all cases regardless 
of whether either mergers or a bar were included. (The slight upturn in the "with mergers" 
line in the plot was not reproducable and is likely a numeric artifact.) Figure [5.191 gives the 
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Figure 5.17: Central relaxation time for the lowest-mass (IMq) stars in the 7 = sphere 
models with a Kroupa IMF. All have A ~ 0.05. 

full density run p{r) at various times in the simulation; interestingly the stellar density rises 
with time at both small and large r, and drops slightly with time in between; this compares 
to the Plummer sphere behavior in which the stellar density dropped with time for the 
largest values of r, e.g., as in Fig. 15.61 The difference could be due to the 7 = sphere 
having relatively more mass to absorb angular momentum from the bar at intermediate 
values of r, to the slower pattern speed of the bar ~ 160 Myr^^, compared to il;, ~ 400 
Myr~^ for model E2A), or simply to the short evolution time of the simulation not allowing 
for angular momemtum to be transported all the way to the cluster's edge. 

As was the case for all the Plummer sphere models, the calculated rate of 25OM0 star 
production remained nill by the end of the simulation, and no stars of mass 125Mq had 
been created via collisional mergers. Table 15.21 shows that model G2A exhibited reverse 
mass segregation (5 < 1); a test run with /max lowered from 3 to 2 showed that S initially 
decreased below 1 and then increased again at times later than were reached by the "full" 
^max = 4 simulation, indicating that the apparent reverse mass segregation is a temporary 
phenomenon. Again, this behavior was also as was found for model E2A. 

5.3.2 Models G3A (Larger Galactic Spheroid) and G3C, Kroupa IMF 

Model G3A started with a somewhat higher density than model G2A did, and had a cor- 
responding faster and greater amount of evolution of the system; model G3C continued 
this trend. In Figures 15.201 and 15.221 the runs without a stellar bar were either unstable or 
ended early, but the case with a bar case continued long enough to reach the beginnings 
of core collapse; these were the only Kroupa-IMF models studied to do so. Even so, the 
density of high-mass stars remained low enough, and the total simulated time was short 
enough, that no 125M0 stars were created by the simulations' end. (The calculated value 
was ~ 0.6 new stars of mass 125M0 for both models, and the Fokker-Planck code ignores 
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Figure 5.18: Central density v. time for 7 = sphere model G2A with a Kroupa IMF, 
^max = 4 and A = 0.051. The nonrotating (A = 0) case gave similar results but was less 
stable and not plotted. 




Figure 5.19: Density p{r) v. radial distance r from the cluster center for the start, midpoint 
and end of the simulation for model G2A with a Kroupa IMF. Neither rotation nor stellar 
mergers were enabled, and /max = 4; this corresponds to the solid-line plot of Fig. 15.181 
The behavior of p{0) v. r for the other two cases from Fig. 15.181 is almost identical to that 
plotted here and are so not shown. 
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Figure 5.20: Central density v. time for 7 = sphere model G3A with a Kroupa IMF, 
^max = 3, A = 0.046 and a 7.7% stellar bar. The gap in the no-merger plot is due to a glitch 
in the output routine. Runs without a stellar bar did not yield stable simulations and are 
not plotted. 

stellar population numbers that are less than 1.) Similarly the calculated rate of 25OM0 
star production remained nill, as given in Table [5^21 A modest amount of mass segregation 
did take place even in the short simulation time. 

Despite the shorter simulation time, the run of p{r) given in Fig. 15.231 displays an increase 
in density for model G3C at small and large r values (and a decrease for intermediate r) 
more clearly than does Fig. 15.191 for model G2A or Fig. 15.211 for G3A. This is consistent 
with G3C's much-shorter relaxation time, e.g., as shown in Fig. 15.171 However, while model 
G3A may be seen as representing a larger version of the galactic spheroid modeled by G2A, 
it is doubtful that model G3C corresponds to a realistic astronomical system, and is more 
useful here as a demonstration case. 

5.3.3 Model G2A, Arches IMF 

Model G2A with an Arches-style IMF behaved much like model G2A/Kroupa but with a 
more rapid evolution. Figure [5.241 shows the central density p(0) increasing roughly linearly 
with time but only a factor of ~ 8 before the simulation ended due to numerical difficulties. 
The presence or absence of either a stellar bar or collisional mergers again had very little 
effect on the overall evolution of the system. Table 15.31 shows that the system starts with 
just a high enough rate of collisions of massive stars G250 = l.lMyr^^ to produce a small 
number of 25OM0 stars, with the rate increasing modestly bebfore the simulation ends. 
And in the short simulation time exactly one star of mass 125M0 was created through 
collisional mergers of smaller stars. 

This chapter has presented the simulation results for each model studied in turn. In the 
next chapter these results will be discussed with respect to what astrophysical effects were 
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Figure 5.21: Density p{r) v. radial distance r from the cluster center for the start, midpoint 
and end of the simulation for model G3A with a Kroupa IMF, A = 0.046, a 7.7% stellar 
bar, and stellar mergers enabled. This corresponds to the dashed plot of Fig. 15.201 




Figure 5.22: Central density v. time for 7 = sphere model G3C with a Kroupa IMF, 
^max = 2, and A = 0.046. The nonrotating (A = 0) case did not yield a stable simulation 
and is not plotted. 
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Figure 5.23: Density p{r) v. radial distance r from the cluster center for the start, midpoint 
and end of the simulation for model G3C with a Kroupa IMF, A = 0.046, a 7.7% stellar 
bar, and stellar mergers enabled. This corresponds to the short-dashed plot of Fig. 15.221 
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Figure 5.24: Central density v. time for 7 = sphere model G2A with an Arches-style 
IMF, /max = 3 and A = 0.051. Two of the runs continued beyond t = 0.07 Myr but became 
unstable; those portions are not shown. 
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observed for which models, how that compares to what would be expected from timescale 
and density arguments, and what conclusions can be reached regarding the production of 
massive central objects in actual stellar systems. 
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Chapter 6 



Discussion 



According to Boker, nuclear clusters - massive, dense stellar clusters in the nuclei of galaxies 
- are observed to be nearly ubiquitous and share a similar relation to the host galaxy as do 
active galactic nuclei [TT|. He closes with the following quote: 

It has recently been proposed by [Ferrearse et ah 2006] that nuclear clusters 
extend the well-known scaling relation between the mass of a galaxy and that 
of its central super-massive black hole (SMBH) to lower masses. This has trig- 
gered speculation about a common formation mechanism of nuclear clusters and 
SMBHs, being governed mostly by the mass of the host galaxy. The idea put 
forward is that nuclear clusters and SMBHs are two incarnations of a central 
massive object which forms in every galaxy. In galaxies above a certain mass 
threshold (~ 10^*^ Mq), galaxies form predominantly SMBHs while lower mass 
galaxies form nuclear clusters. 

... Is a nuclear cluster possibly a pre-requisite for the formation of a SMBH? Is 
the formation of a BH (not necessarily a super-massive one) a logical consequence 
of the high stellar densities present in nuclear clusters? Progress along these lines 
will require a better understanding of the formation of pure disk galaxies in the 
early universe, as well as improved models for the evolution of extremely dense 
stellar systems. 

This study is an attempt at just such an improved model of dense stellar systems, linking 
what the above reference calls nuclear clusters with the eventual supermassive black holes 
in the centers of active galaxies. This chapter begins with two brief discussions, first of 
the timescale on which mass segregation alone would be expected to result in core collapse, 
and then of what central density is required before either runaway accretion or runaway 
collisional mergers might occur. The main body of the chapter then discusses the physical 
effects exhibited by the various simulations and what they imply for astronomical systems. 
Finally a brief closing considers what future observations may support the assumptions and 
results of these simulations, and what further improvements of the simulation method are 
indicated. 



96 



6.1 Timescale Arguments and Expectations 



6.1.1 Mass Segregation 

For mass segregation to occur, energy equipartition across different mass populations must 
be an unstable situation. Spitzer in 1969 formulated an analytic condition for preventing 
the mass-segregation instability: for a system of total mass Mi + M2, with Mi contained 
in stars of individual mass mi and M2 in stars of individual mass m2 > mi as 



which for realistic mass spectra is never satisfied, i.e., there are enough heavy stars that 
they decouple from the dynamics of the system as a whole and experience mass segregation 
|54j . For a system with multiple masses, the analogous condition has been numerically 
determined to be [H] 



in which iV; is now the total number of stars of individual mass mi and is the total 
number mass mh > mi. For the Kroupa IMF used here, the value of the above relation 
is 13; for the Arches IMF, 209. Even if the Kroupa IMF is extended down to 0.2Mq, the 
value is still 1.4. So indeed for all IMF studied here there are sufficient heavy stars so that 
mass segregation should at least be possible. 

Given the plausibility of mass segregation, one can ask if it would be expected to foster core 
collapse. The timescale for segregation-driven core collapse has been variously estimated 
to be ~ 10% of the initial half-mass relaxation time [41] or ~ 15% of the initial mean-mass 
central relaxation time [54j, both independent of the IMF used. For the models studied only 
the extreme model G3C has a resulting expected segregation-collapse timescale which falls 
under the 3 Myr stellar- lifetime limit, being in the 1 — 2 Myr range for the two conditions 
given above. (Figure [5. 171 shows model G3C having a very short relaxation time in general.) 
The runs of model G3C all ended before 0.5 Myr of simulated time and so would not be 
expected to exhibit segregation-driven collapse. In addition, other models were allowed to 
run beyond the 3 Myr limit as a test but none reached its expected segregation-collapse 
time - the closest being model E4B with an Arches IMF and no stellar bar, which had an 
expected collapse time of tc — 15 Myr and which ended at t ~ 10 Myr showing no signs of 
collapse. 

6.1.2 Critical Density 

Whether or not core collapse is achieved, the central question is whether, given the ~ I5OM0 
upper mass limit of the IMF, a more-massive object can be formed which will evolve to be 
a seed black hole. Once a critical local density of p > pcrit = 5 x 10^(yq(^^^)Mq/pc^ is 
achieved, a ~ IOOMq star can directly accrete W^Mq of material in ~ 5 Myr, creating an 
IMBH [55J. As all models used here have a central velocity dispersion cji > 100 km/s, Tables 
15.21 and 15.31 show that none of the models meet that requirement. Model G3A/Kroupa may 
eventually reach it once deeper into core collapse than the simulation was able to track: Fig. 




(6.1) 




(6.2) 



97 



15.161 shows its ai ~ 150 km/s and remaining roughly steadly even as core cohapse starts 
at t ~ 0.4 Myr in Fig. 15.201 Model E4B/ Arches reaches an even higher central density 
during its core collapse (Fig. 15. 4p . but it has a central velocity dispersion approaching 500 
km/s and still increasing (Fig. 15. Sp . However neither model is expected to represent a 
typical astronomical system. As described in ^2.61 G3A has a somewhat high density for 
the galactic spheroids that the 7 = spheres best model. Likewise the Arches IMF used 
in model E4B/Arches is determined from stellar clusters, while E4B's velocity dispersion 
places it in the range of galaxy cores, not individual clusters. 

Other than accretion, collisional mergers are another avenue for forming massive objects 
from smaller stars. For massive {12bMQ) stars, a local density of p^^^ ~ 1.2 x 10^Mq/pc^ is 
required in order for the average star to experience I collision per main sequence lifetime |88j . 
None of the models started with such a high central density, and only model E4B / Arches 
with the maximally strong bar achieved it by the end of the simulation. While by no means 
is it necessary for a collision rate of 1 per 125M0 star be reached in order for a massive 
object to form - such a rate would produce an abundance of 250Mo objects, while previous 
simulations of young clusters with an upper mass limit of I2OM0 have found that more 
than one very massive star is never formed [56j - it does give an indication of the density 
range required. As will be seen in the following discussion of detailed simulation results, for 
a broad range of initial models using the Arches IMF (observed in stellar clusters), 25OM0 
objects are calculated to have been formed, while with the Kroupa IMF (representing the 
field star population) no such massive objects are created. 

6.2 Summary of Simulation Results 

The models studied fall into two main classes. The majority are Plummer-sphere initial 
models, which to be consistent with Quinlan and Shapiro [7j are labeled with an initial "E". 
Models E2A and E2B can represent either dense globular clusters or the nuclei of dwarf 
elliptical or bulgeless spiral galaxies, and E4B the core of a giant elliptical galaxy. 

The other class of model is the "7 = sphere", labeled with a "G", which fit the surface 
brightness of galactic spheroids. Model G2A has identical mass and core radius as model 
E2A, while models G3A and G3C have the same mass as model E2B but with smaller 
core radii. G3C is an extreme case which does not likely represent a realistic astronomical 
system. 

For each model two possible initial mass functions were employed: a Kroupa IMF based on 
general field star populations, and an Arches-style IMF similar to that observed in dense 
clusters specifically. Note that the Arches IMF may not be a realistic IMF for either model 
E4B or the 7 = spheres, as those models best represent galactic centers while the Arches 
IMF is derived from observations of individual stellar clusters. However, the simulations 
of those models serve to complete the overall picture of the physical evolution of stellar 
systems at galactic centers. More complete details of the individual models and IMFs, 
and how rotation is introduced into them, is given in ^2.61 The remainder of this section 
describes the physical effects exhibited by the various simulations, first for the Plummer 
models and then for the 7 = spheres. 
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6.2.1 Plummer-Sphere Simulations 



The simulations which started with a Plummer Sphere distribution (the "E" models) give 
a clear picture. As each model's results show, the presence of a moderate-strength stellar 
bar can dramatically raise the rate of central density increase of a dense stellar cluster, by 
a factor of 1-2 orders of magnitude. However, the IMF also plays a crucial role: with the 
standard Kroupa IMF the rate of increase of central density p{0) remains roughly linear 
with time through the main-sequence lifetime of the largest stars in a cluster, even for the 
most-dense model studied (E4B) - as shown in Fig. 15.11 - whereas with the cluster-specific 
Arches-style IMF core collapse can take place near to the 3 Myr stellar lifetime limit {e.g., 
model E2B with a 5% bar, shown in Fig. 15. lip , or even within it (model E4B with an 8 or 
10% bar. Fig. EM- 

More detailed discussion of the various physical aspects is given below. 
Rotating v. Non-rotating Clusters 

As described in ^6.1.11 none of the Plummer models was expected to experience mass- 
segregation-driven core collapse, and none was observed in the simulations. For all of 
models E4B/Kroupa (Fig. [53D, E4B/ Arches (Fig. [53D and E2B/ Arches (Fig. [5lT|) both 
the nonrotating and rotating cases showed only slow growth in the central density by sim- 
ulation's end, if no stellar bar was incorporated. Interestingly, rotation was not a strong 
impediment to this essentially secular increase in p(0), e.g., the E4B/Kroupa cases' density 
plots track each other with the only difference being the lower initial p{0) value of the rotat- 
ing cluster, while in the E4B/ Arches and E2B/ Arches models the rotating case actually has 
a somewhat larger rate of central density increase. While counterintuitive, rotation serving 
to enhance the rate of collapse has been observed before in Fokker-Planck simulations |39j, 
possibly due to individual stars being able to interact with a larger (or at least different) 
range of neighboring stellar orbits. 

Effects of a Stellar Bar Perturbation 

When a moderate stellar bar ~ 5% to 10% strength, in line with observed bars - was 
included, however, all models showed rapid increase in central density. Core collapse was 
achieved in both E4B/Arches and E2B/Arches, and in less than the 3 Myr stellar-lifetime 
limit when the bar strength was ~ 8% or greater. Figure 15.131 shows the dramatic effect a 
greater bar strength can have on reducing the time to core collapse, going from t ~ 4 Myr 
with a 5% bar to slightly over 2 Myr with an 8% bar. Comparing Figures 15.101 and 15.111 it 
appears that bar-driven core collapse would also be expected in model E2A/ Arches had its 
simulation been able to remain numerically stable for a longer amount of simulated time. 

Collisional Stellar Mergers: Effect on Cluster Dynamics 

In general stellar mergers had little effect on the overall dynamics in the Plummer-sphere 
simulations. Figures l5.1l for E4B/Kroupa and 15.111 for E2B/ Arches show only a very slight 
extra increase in p{0) when stellar mergers are enabled, as is also born out in comparing 
values of S and S' in Tables 15.21 and 15.31 Even in the most extreme Plummer-sphere model 
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studied, E4B/ Arches with a 10% stellar bar, the main effect was that core collapse was 
achieved slightly earlier with stellar mergers, as shown in Fig. 15.41 That figure does hint 
that as core collapse progresses, collisional mergers may begin to drive the evolution of the 
central regions. However, this computational model cannot track into that regime, and it 
remains that the core collapse that is seen here is consistently bar-driven. 

Collisional Stellar Mergers: Massive Star Formation 

When the buildup of massive {250Mq) stars via collisional stellar mergers is considered, the 
accelerated central density increase brought about by the bar perturbation is insufficient to 
produce even a single massive star within 3 Myr using the Kroupa IMF, as seen in Table [5^ 
But the Arches IMF is already sufficiently populated with high-mass (up to 125Mq) stars 
to allow for the collisional production of massive 250Mq stars, and the accelerated collapse 
induced by the bar only increases the rate of production, as seen in the 6*250 (ii) column 
in Table 15.31 So as far as the formation of massive stars is concerned it is the IMF which 
determines whether or not 25OM0 or larger objects are created at all. The overall rate of 
formation due to collisional mergers is primarily set by the initial conditions of the cluster, 
i.e., its size and density, as can be seen by comparing the ^250(0) values for models E2A, 
E2B and E4B in Table 15.31 However, the table also shows that the presence and strength 
of any stellar bar can have a large effect on how much the rate of massive star formation 
G250 increases in the period before core collapse occurs. 

Complementing the above discussion of the calculated rate of massive object (250Mq) 
formation, one can look at the numbers for the largest stars in the IMFs used. The physical 
upper limit was I5OM0, which was implemented numerically as a mass bin centered on and 
associated with stars of mass 125M0. In the Kroupa IMF models, no simulation yielded 
an increase in the number of 125M0 stars due to collisional mergers, which is consistent 
with a null rate for 25OM0 object formation. As detailed in Chapter [5l the Arches-IMF 
models ranged from creating 15 stars of mass 125M0 stars in model E2A with a 5.3% bar, 
to 500 stars of mass 125M0 stars in model E4B with a 10% bar. Noting the short amount 
of simulation time for which model E2A ran, these numbers are commensurate with the 
calculated rates of 25OM0 object formation. 

Collisional Stellar Mergers: Possibility of Runaway Merging 

The question of how massive-star production due to collisional mergers might be self- 
reinforcing is more difficult to answer. It would be expected that as more low-mass stars 
merge and so produce more intermediate-mass products, then more intermediate-mass 
merger candidates would be available to again collide and produce even more-massive stars. 
As discussed in ^5.2.41 model E2B/ Arches was allowed to run until time t = 3 Myr without 
stellar collisions enabled, and the rate of massive-star production due to mergers was at 
that point calculated to be within 3% of the value found when mergers had been allowed 
ab initio (and both were > 25% larger than the rate at t = 0). Thus the amplification of 
the rate of massive star production due to the collisional merging of smaller stars is not 
significant when compared to the increase in the merger rate which results from the overall 
central density increase afforded by the stellar bar, at least in the pre-core-coUapse period. 
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6.2.2 "7 = Sphere" Simulations 

The 7 = sphere models have similar density profiles to their Plummer sphere counterparts, 
although they are somewhat less centrally concentrated. A major distinction is that the 
7 = spheres have a much shorter initial relaxation time, as seen by comparing Figures 
15.171 and 15.91 This results in a very different physical evolution, as described next. 

Rotating v. Non-rotating Clusters, Effect of Bar Perturbations 

Contrasting the situation for the Plummer-sphere clusters, simulations of initial 7 = 
spheres showed no significant differences in evolution between nonrotating systems, rotating 
systems without a stellar bar, and rotating systems with a bar. Figures 15.181 for model 
G2A/Kroupa and 1 5. 241 for G2A/Arches show this with respect to the central density p(0). In 
model G2A/Kroupa the with-bar case exhibits a somewhat more rapid increase in p{0), but 
it is extremely slight compared to the effect the bar has in the Plummer sphere simulations. 
In Fig. 15.241 model G2 A/ Arches shows even less difference, although it ran for a very short 
simulated time. 

Even though the simulations of 7 = sphere models were not able to track the system as 
far into their evolution as some of the Plummer-sphere simulations were able to do. Figures 
15.181 (for G2A/Kroupa) and 15.221 (for G3C/Kroupa) show a trend for the 7 = sphere 
models with no stellar bar to evolve almost identically to the same models with a bar, in 
stark contrast to how the Plummer sphere models behaved. Thus for the 7 = sphere 
models, overall system relaxation dominates over any bar- or collision-driven collapse, and 
again rotation does not obviously support the system against an increase in the central 
density. 

CoUisional Stellar Mergers 

Stellar mergers played a similar role in 7 = sphere simulations as they did for the Plummer- 
sphere models. Figures 15.181 and 15.201 show that the central density p{0), in cases in which 
stellar mergers are enabled, exactly tracks p{0) in the corresponding no-mergers simulations 
for models G2A and G3A with a Kroupa IMF. Mergers have only a slight effect on p(0) 
for model G2A/Arches (Fig. I5.24p . Only in the extreme model G3C of Fig. 15.221 do stellar 
mergers affect the evolution, and then only by producing a slightly more rapid start to core 
collapse. 

In terms of numbers of high-mass stars formed through collisions, the 7 = situation is 
slightly different, in that for both models G3A/Kroupa and G3C/Kroupa the calculated 
number of 125M0 stars formed via collisions of lower-mass stars was between 0.5 and 1. 
Given the short simulation time of each model, a few additional 125M0 would be expected 
to be created within the 3 Myr limit. While not interesting astrophysically, this contrasts 
with the situation for Plummer spheres in which no Kroupa IMF simulation produced any 
125Mq stars. However, Table [5^ still shows that no 250Mq objects at all would be created 
through mergers even in these cases - and Table 15.31 shows that only a small number are 
formed even in model G2A with a more top-heavy Arches IMF. 

What remains unclear is whether model G3A/Kroupa will reach a regime in which massive 
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star formation becomes plausible: as seen in Fig. 15.201 it does achieve the beginnings of 
core collapse, unlike any Plummer model with a Kroupa IMF. It does so in a short amount 
of simulation time, leaving ~ 2.5 Myr for post-core collapse evolution to occur before the 
stellar population begins to evolve off the main sequence. Model G3C/Kroupa achieves 
core collapse even more rapidly, but it is doubtful that it represents a realistic astronomical 
system. 

6.2.3 Summary 

Before performing the simulations it was expected that rotational support against collapse 
would be countered by the outward transport of angular momentum afforded by a stellar 
bar, allowing core collapse to proceed even in a rotating system. While the effect of the 
bar was observed in all Plummer sphere simulations it did not result in the anticipated core 
collapse in all models, nor in an ensuing runaway of collisional mergers creating a massive 
object. For the 7 = spheres the overall system relaxation dominated over the effect of the 
bar. 

In the end the initial mass function turned out to be the determining factor: with the 
possible exception of model G3A no simulation performed using a Kroupa IMF, represen- 
tative of a general galactic stellar population, reached core collapse, and none obtained a 
large enough stellar merger rate to create a massive object. In contrast, all models which 
employed the Arches-style IMF, which is specific to dense stellar clusters, started with a 
merger rate already sufficient to produce massive 25OM0 objects - although the presence 
of a stellar bar did make the difference in whether or not core collapse was reached, as well 
as in producing a large increase in the massive-object formation rate. 

The above description is fairly robust: 

• with a Kroupa IMF only the most extreme models (G3A and G3C) reached core 
collapse; 

• with an Arches IMF all models reached core collapse (or for E2A, showed indications 
of it); 

• no 7 = sphere model was dominated by the effect of a bar; 

• even with an Arches IMF, a bar was required in order for any Plummer model to 
reach core collapse; and 

• runaway mergers were not observed in any simulation, although a large increase in 
high-mass merger rates followed from core collapse in barred models. 

Of note is that a cluster's evolution did not change qualitatively even when the merger rate 
was artificially increased by a factor of 5, as described in §5.2.61 - so it does not appear 
that a collisional-merger runaway is likely to occur for any pre-core-collapse system of main 
sequence stars. 

Thus the results indicate two possible paths for formation of a massive object within a dense 
stellar system, one clear and one somewhat tentative. For cluster-sized systems similar to 
models E2A or E2B - including dwarf elliptical galaxy nuclei or bulgeless spiral galaxies - 
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with an Arches-style IMF, the expected degree of rotation is sufficient to support a stellar bar 
perturbation which in turn leads to core collapse through the outward transport of angular 
momentum. However, even without core collapse, through stellar mergers the system will 
produce massive objects in less than the main-sequence stellar lifetime sufficient to seed 
supermassive black hole growth. For larger systems, e.g., galactic spheroids similar to model 
G3A, core collapse can be obtained even with a Kroupa IMF, although massive object 
formation through stellar mergers is not expected to be achieved within main-sequence 
lifetimes and so will require interactions of collapsed relativistic objects, which was outside 
the scope of this study. 

6.3 Future Observations, Future Work 

A basic assumption in this study is that stars preceded quasar black holes in the early 
universe. Two future space missions currently being considered by the European Space 
Agency in its Cosmic Vision 2015-2025 program will go a long way towards addressing this 
assumption by observing the earliest luminous objects whose surrounding clouds of dense 
gas and dust obscure them from current instruments. First, the Far Infrared Interferometer 
(FIRI) will perform high-resolution imaging spectroscopy to resolve the creation of the first 
luminous objects in the universe, separating out the formation of stars and the growth of 
black holes [89]. Complementing FIRI, the ESA's XEUS X-ray observatory will detect the 
earliest quasars as they form. Together these two missions, or similar ones, could answer the 
question of whether stars or quasar black holes formed first [90] . They or similar observations 
may also address the main question remaining about the current study's results: whether 
or not an Arches-style IMF, with its weighting towards higher stellar masses, is in fact a 
better representation of the actual IMF that existed in early-universe stellar clusters than 
the Kroupa IMF is. If so, then the main result of these simulations - that individual nuclear 
stellar clusters are likely to form seed objects for massive black holes - holds. If not, then 
one must resort to post-core-collapse evolution to obtain the seed objects, and likely in 
galactic spheroids instead of in nuclear clusters. 

As with all simulations, the results presented here have some limitations that could be 
addressed in future work. 

• Finite stellar lifetimes. A conceptually straightforward but technically nontrivial 
extension would be to track the system's dynamical behavior as the more-massive 
stars move off the main sequence, which would involve the appropriate moving of 
stars from their present distribution functions to new ones representing white dwarfs, 
neutron stars and stellar-mass black holes - as represented by terms Bg and Rg in the 
general Fokker-Planck equation of 13.51 but not yet implemented in the simulations. 
Corresponding different collisional cross sections would also be required, as has been 
done by Quinlan and Shapiro in modeling one-dimensional, nonrotating Plummer 
sphere systems ([7J, [Ij). 

• Rotation-induced flattening. As described in §2.51 the provision for a non-spherical 
gravitational potential is already implemented in the simulation code, but it is not 
yet debugged and shown to be numerically stable; once that is achieved it may be 
possible to add another degree of realism to the simulations of rotating systems. For 



103 



globular clusters at least, flattening of the system has been found to correlate with 
rotation [251. and to increase merger rates somewhat [26j . which could strengthen the 
results presented here. 

• Dark Matter. No provision has been made for there to be a background component 
in the gravitational potential, e.g., due to the presense of collionless dark matter. This 
is consistent with all previous Fokker-Planck studies of stellar systems - in recent 
works Fiestas et al. ([25j), Kim et al. ([39], [24]) and Takahashi ([45], [46]) all 
ignore dark matter in their globular cluster simulations, while Arabadjis [26] and 
Quinlan & Shapiro (|lj, [7j) look to the collisionally-produced massive objects in their 
simulations as being dark matter. Amaro-Seoane [H] points out that the current 
state of astrophysical modeling "comes to its limits" in embedding simulated systems 
in collisionless dark halos. Still, it is conceivable that introducing an ad hoc offset 
to the overall gravitational potential, or an additional collisionless component to the 
system's matter-density profile, could serve as a toy model of the effect of a smooth 
dark matter background. 

The early evolution of quasars and host galaxies is not simple ~ we have only observed 
the final ~ 30% of quasar evolution, but even that is enough to indicate galactic spheroids 
and massive black holes do not grow together pO]. The results presented here show that 
even when rotation is accounted for, dense clusters of main-sequence stars in early galaxy 
cores can provide a possible mechanism for forming the seed objects that eventually become 
quasar black holes. 
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Appendix A 

Non-Spherical Gravitational 
Potentials 



A.l Potential Calculations in Homeoidal Coordinates 

For the more general, but rarely needed case in which there is an overall and variable 
ellipticity e(a) to the potential, a somewhat more complex procedure than that given in 
^2.51 is required. One still starts with the inverted Poisson equation as given in (j2.40p but 
now the derivation of the inverted Laplacian is more involved, as developed below. 

To describe the non-spherical potential ^ when the isodensity surfaces of p are ellipsoids, 
homeoidal coordinates are called for [42j. Consider the potential of a thin ellipsoid of 
mass dM, semimajor axis a and ellipticity e at homeoidal radius Uq- Defining A = ae, the 
homeoidal radius u is related to the standard cylindrical coordinates by 



cosh^ u sinh^ u 

(There is also an angle-like coordinate v which does not come into play.) Then the potential 
d^{u) at u due to the ellipsoid at Uo is 

\ GdM ( sin^^e, u < Uq , k ^\ 

d^{u) = — <^ ■ -If . ^ ^ (A.2) 

A [ sm (sechnj, u>Uo 

where the mass element dM = pdV for the ellipsoidal shel0 is 

pdV = 47r/ja2(l - 6^)^/2 (A.3) 

and the total potential is found by integrating over mass shells at all values of Uo as will 
now be shown. Let us define d^e = —Q^L ^{^-^ g and d^^u = — ^^^^ sin~"'^(sech n). Having 



^Here it is implictly assumed that the ellipticity e is a slowly-enough varying function of radial coordinate 
a so as to not affect the differential dV; this approximation is examined in the Appendix. 
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no explicit dependence on u, the integration of is straightforward: 



$e(n) = -47rG H (1 - 6^)1/2^ ^ia. (A.4) 

Ja{u) oe 



For the contribution to <I>(n) due to shehs with Uo < u we need a relation for u{a). Noting 
that in the z = plane, i? = a by definition, we accomplish this by solving (jA.ip for u; the 
resulting quadratic equation in sech^ u yields 



sech^ u 



1 



where the necessity for sech^ u < 1 implies that only the upper sign is allowed {e.g., taking 
R = z makes this clear). Then the contribution to the total potential due to mass shells 
interior to u is obtained by integrating ()A.2p using (|A.3P : 

<^uin) = -4nG - e^f- p'^^^^^^ da. (A.6) 

Jo ae 

Equation (jA.SP may be expanded, yielding approximations necessary for use when R — t- 
0. (When A — t- the simpler spherically-symmetric procedure can be used.) The total 
gravitational potential at an arbitrary value point u{R, z) is thus given by 

(D-w^^) = + (A.7) 

Note that with use of (jA.SP in ()A.6p . the only reference to the homeoidal coordinate u in 
the calculation of <I>(n) is in locating the transition from to <!>£) despite Uq being the 
(implicit) variable of integration. Note also that the coordinate transformation of (jA.ip is 
referenced to the integration variable Uq, and not to the location of measurement u; thus 
each homeoidal transform is defined by the mass shell responsible for the element of 
currently contributing to = J d^u- Hence a = a[uo) and e = e{uo)] neither a nor e 
has an explicit dependence on u. Likewise, the A in (lA.Sp is A = a{uo)e{uo) even though 
the coordinate being transformed is u and not Mq, and so (jA.ip constitutes a "running 
coordinate transformation" which is continuously evolving across the range of integration 

of im . 



A. 2 Laplacian Solution in Ellipsoidal Coordinates: Numeri- 
cal Aspects 

A novel aspect of this study is that e is allowed to vary over the range of a. This is different 
from the fixed-ellipiticy methods {e.g., as described by Tremaine and Weinberg [30] ) but was 
found to be both advantageous in that it allows potentially more accurate modeling, and 
necessary to avoid a pitfall of the fixed-ellipiticy schemes: physically interesting interactions 
will primarily occur in high-density regions, i.e., the core of the distribution, but these tend 
not to be the high-ellipticity regions which are usually far from the center of the cluster 
~ and so taking a single value of e for the entire cluster can artifically increase the effects 
of non-sphericity. This was found to be the case in practice, and was what motivated the 
allowance for a variable e. (Provision for fixing e at a predetermined value is also built in 
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however.) 

One may worry that a changing e(a) may induce troublesome shell-crossing. Numericahy, 
this is only a concern if the a} grid is so fine and e(a) changes so rapidly that the semiminor 
axes h = a(l — e)^/^ of adjacent shells on the grid cross (i.e., if ar+i(l — &r+\f'^'^ < ar(l — 
e^)^/^ for gridpoints and a^^^ > a^). In practice the variability of e(a) is not large 
enough for this to occur, given the courseness of the grids. And in principle, possible 
shell-crossing is not a fundamental concern anyway, as the ellipsoidal isodensity surfaces are 
merely smoothed representations of the underlying discrete stellar distribution. 

To see how much effect the variability of e(o) has on the entire p — $ scheme, consider what 
effect it has on the mass element dM, which was not taken into account in ()A.3p . This is 
not a strict calculation of what effect variable e has on the calculus, but it does give an idea 
of the size of its effect: 



dM = pdV = pd 



A-Kpa^l - e^) 



2^,1/2 



ae de 
^~ 3(l-e2)i/2d^ 



da (A. 



which shows that the effect of the dependence of e on a is of 0{e'^) smaller than the dom- 
inant term. Using the more straightforward form of (jA.Gh greatly simplifies the numerical 
calculation without explicitly affecting the results, as the e(a) run is still allowed to converge 
on a consistent set of values along with p{a) and ^>(a) during the iterative solution for all 
of these at each new timestep, with resulting new /. 

Thus the assumption of ellipsoidal isodensity surfaces, along with /„, determines the poten- 
tial The required integrals described above, however, are too computationally expensive 
to be performed each time the knowledge of potential is needed in the Fokker-Planck coeffi- 
cient calculation. After testing interpolation schemes and analytic approximations, only the 
Glutton-Brock self-consistent field ( "SCF" ) method as described by Hernquist and Ostriker 
[91] proved adequate in both accuracy and speed. For the SCF "core radius" rcore the value 
at which $(rcore) = *&(0)/\/2 is used. This is also the core radius of the Plummer potential 
with central potential $(0). The field method is tested at each timestep, and if sufficient ac- 
curacy cannot be achieved with a reasonable number of expansion terms, the code falls back 
on direct integration. For interpolating over the other, one-dimensional grids of quantities 
discussed in Chapter El simple polynomial or cubic spline schemes are employed. 
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Appendix B 

Tables of Symbols 
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Quantity Symbol Example 

^2.21 radial component r or 1 Vr, Ii 

tangential component T or 2 vt, I2 

azimuthal compoment 2; or 3 J^, 

orbital endpoints p and a rp, TZa 

dummy component variable i oi j = 1, 2, or 3 ilj, Dij 

^2.2.21 "circular" value or component c Tc 

§2.3.21 maximum dynamically- allowed value max Jmax, /3max 

§2.3.31 velocity-space coordinates v v, 6y, cp^ 

a specific mass "bin" q; later also q' k, Q niq, fg' 

§2.4.11 total value of a dynamical quantity tot Jtot 

minimum dynamically- allowed value min ^min 

§2.4.41 average over a given quantity x, 

possibly with a weighting function (•)^ see text for definitions 

§2.61 core (radius) core rcore 

^2.6.21 amount of a dynamical quantity 

due to rotation rot Jrot, Trot 

^3.2.11 spherical harmonic indices m 

resonance indices knm; later ^i^2^3 ^/cnm 
an individual object (whether 

a single star or the stellar bar) * 

§3.31 the stellar bar B fs, 



^3.41 repeated summation index [any repeated index] £pO,p = ip^p 

Table B.l: List of frequently- used subscripts. 
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Quantity 



Symbol 



mil 
mil 



Newton's constant of gravitation 

radial component of orbital action 

tangential component of orbital action 

vertical component of orbital action 

orbital energy per unit mass 

orbital position coordinates 

radial &: tangential orbital frequencies 

gravitational potential per unit mass 

canonical orbital angle coordinates 

orbital inclination (Euler angle) 

orbital azimuth (Euler angle) 

orbital radial endpoints 

radial component of orbital velocity 

"circular" orbital radius 

individual stellar mass 

tangential component of orbital velocity 

average rotational velocity 

orbital polar angle 

"ellipsoidal" radial coordinate 

cylindrical radial and vertical coordinates 

one-dimensional (nonrotating) 

distribution function for stars of mass niq 
total cluster gravitational potential energy 
total cluster rotational angular momentum 

cosmological rotation parameter 



G 

11 or /; later also I 

12 or J or J; later also J' 
Jz = J cos j3 

E 

r , 0, 4> 
<^ 

wi, W2, If 3; later also Wi etc. 

/3 



' pi 



ra, later also TZp, TZa 
vr = V[2(^-$)- j2/r2] 
Tc where dvr{rc)/dr = 
m 

vt = J/r 
w{r) 

a = f -/3 
a(r, e(r)) 
R, z 



J, 



grav 
rot 



A 



Table B.2: Symbols first introduced in Chapter [21 listed by section in which each was first 
introduced. These denote properties of individual stars in the system, not bulk properties. 
Bulk properties based on summing over the individual stars' values are denoted by sub- 
scripts, for example total energy of the stellar system E^ot or the rotational contribution 
to the total angular momentum Jfot- Only symbols that are used in more than one section 
are included. 
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Quantity 



Symbol 



[2:3:2] 

EXH 



gravitational potential per unit mass $ 
stellar mass density p 
total cluster mass M 
stellar distribution function in position space F(r, 
stellar distribution function 

in energy/angular momentum space 
orbit-averaged stellar distribution function 
stellar distribution function parameterized 

in terms of orbital inclination 
even {g) and odd (0) components 

of inclination parametization 
total angular momentum of all stars of energy E 
ellipticity of stellar mass distribution 
one-dimensional "rms" stellar velocity dispersion cji 
anisotropy parameter 5 
generalized stellar velocity dispersion ao = \/(l 



F{E,J) 

f^{I,J;(3)^h{amf{I,J) 

h{a) = g{\a\) + e{a) 
JtotiE) 
e(r) 



6)al 



Table B.3: Symbols first introduced in Chapter [21 listed by section in which each was first 
introduced. These are bulk properties of the system or of a subpopulation thereof. Only 
symbols that are used in more than one section are included. 
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Quantity 



Symbol 



T] Fokker-Planck drift coefficient 
Fokker-Planck diffusion coefficient 
ad hoc merger loss and gain terms 
in F-P equation 
2.1| gravitational potential of individual perturber 
(whether a single star or the stellar bar) 
spherical harmonics 
radial coefficients 

of perturbation potential expansion 
coefficients of perturbation potential 

expansion in action space 
Slater rotation coefficient 
resonance strength coefficient 
resonance indices 
2.2l any quantity y{(3), weighted-averaged 
over orbital inclination 
any quantity y{l3)x{9), weighted-averaged 

over polar angle 9(1^,13) 
average-squared perturbation 

coefficient strength 
perturbation "turning on" parameter 
(complex) perturbation orbital frequency 
2.1 [ Relative density increase of highest-mass stars 
As above, but including effect of stellar mergers 
Rate of massive (25OM0) star creation 



Cj{I,J) 
Dij{I,J) 

Lg{I,J), G,{I,J) 

^ knm (^1 ^) 
VlnmW) 
Wklnm{I, J) 

{k,n,m); later (^i,-^2;-^3) 



{y,x)e = ^Jdl3h{fi)y{P)dijx{e) 



UJ = mQ^ + ir] 

s = [pQ{o)/pmtj[PQ{o)/pmt, 

S' = ... 

G250 



Table B.4: Symbols first introduced in Chapters [3l [4] and [H listed by section in which each 
was first introduced. Only symbols that are used in more than one section are included. 
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